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AN EARTH- MARS MISSION -ANALYSIS PROGRAM 


By James F. iCibler 
Langley Research Center 

SUMMARY 

A rapid, flexible, preliminary Earth-Mars mission-analysis computer program has 
been developed. The program computes a conic interplanetary trajectory approximation, 
a noncoplanar impulsive deboost maneuver into a closed orbit about the target planet, and 
many mission-dependent and mission-independent parameters to allow examination of the 
entire flight profile. The capabilities of the program are discussed along with the require- 
ments for computing a general planet -to-planet mission. Examples of program input and 
output and sample data analyses are presented for an Earth-Mars mission during the 1973 
launch opportunity. A flow diagram for the main program, the input and output descrip- 
tion, and a complete program listing are presented in the appendixes. 

INTRODUCTION 

Interplanetary flights from Earth to the planets represent a significant part of the 
space effort. Detailed study is required for each of these flights. A necessary part of 
the study is a preliminary mission analysis which consists of choosing a suitable mission 
profile from an infinite set of possible candidates. Many variables and trade-offs are 
available to the flight planner. For example, once a rocket booster is chosen, the maxi- 
mum allowable spacecraft weight is specified for a given launch date and arrival date. 

The planner must then trade weight for required launch-arrival periods until a feasible 
combination is obtained. Once the spacecraft weight is determined, an allocation must be 
made for fuel to perform midcourse corrections, a deboost maneuver at the planet, and a 
deorbit maneuver to the surface of the planet. Each of these maneuvers will depend upon 
other mission constraints such as the desired orbit at the planet and the desired landing 
point on the surface. In addition, the mission planner must consider the effect of scien- 
tific requirements on the mission profile. For example, the landing point must be located 
in a scientifically interesting area and must have proper lighting for any onboard optical 
equipment. The orbit about the planet must satisfy constraints such as communication 
requirements with the Earth and the necessity for solar cells to be exposed to sunlight 
for the greater part of each orbit. The many different problems involved in preliminary 
mission analyses present a real task for the flight planner. 



At Langley Research Center, computer programs have been developed to solve sev- 
eral individual parts of the mission-analysis problem. However, experience with the 
Viking project has shown the difficulty of data interchange between the programs and the 
necessity for an integrated approach to a preliminary mission analysis. The program 
described herein is an attempt to combine the many facets of preliminary mission design 
into one rapid and flexible program. In addition, capability not previously available in 
program form, such as a noncoplanar impulsive -burn deboost maneuver, has been included 
in this program. 

The present version of the mission-analysis program is concerned only with Earth- 
Mars missions of the Viking type. However, modifications described herein would allow 
study of interplanetary missions to other planets. The accuracy of the program is lim- 
ited by the use of Keplerian mechanics and impulsive-burn maneuvers rather than finite- 
burn integrating schemes. However, it is felt that for preliminary mission design, the 
order-of-magnitude accuracy involved in the approximations, as compared with an inte- 
grating program, is far outweighed by the several orders of magnitude gained in compu- 
tational speed and program flexibility. Results from the various program elements 
agree with results obtained from other conic programs that were previously designed 
individually to study specific parts of the total mission. 

Information required for operation of the program is contained in the appendixes. 
Appendix A includes a brief flow chart and a description of the primary subroutines. An 
explanation of the required input is given in appendix B. Appendix C describes the output 
parameters, and a complete FORTRAN listing is given in appendix D. The program was 
developed for use on a Control Data 6600 series computer and requires a field length of 
approximately 65000s. 


SYMBOLS 

B a vector from center of target planet and perpendicular to approach 

asymptote of incoming hyperbola 

C 3 twice total geocentric injection energy per unit mass, km^/sec^ 

Dla declination of launch asymptote as measured from Earth’s equator, deg 

f true anomaly, deg 

f^eorbit anomaly of deorbit, deg 

G Sun lighting angle, deg 


ha 

2 


apoapsis altitude of specified elliptical orbit, km 



hp 

i 

J2 

k 

ra 

^entry 

R 

S 

T 

Ve 

Vh 

Voo 

AV 

a. 

M 


periapsis altitude of specified elliptical orbit, km 
inclination of elliptical orbit, deg 
oblateness coefficient for Mars 

true anomaly of landing point (periapsis to landing-point angle on ellipse) , deg 

apoapsis radius of ellipse, km 

radius at entry to Martian atmosphere, km 

periapsis radius of ellipse, km 

radius of Mars, km 

K X T with R a unit vector completing the RST triad 

unit vector, parallel to approach asymptote at Mars and passing through 
center of planet 

unit vector perpendicular to S’ and parallel to ecliptic plane 
velocity on ellipse, km/sec 
velocity on hyperbola, km/sec 

hyperbolic excess velocity of spacecraft relative to Mars, km/sec 

deboost velocity- change requirement, km/sec 

flight -path angle at entry into Martian atmosphere, deg 

declination of landing point, deg 

declination of subsolar point, deg 

gravitational constant for Mars, km^/sec^ 

right ascension of landing point, deg 
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0s 

right ascension of subsolar point, deg 


right ascension of ascending node, deg 

a; 

argument of periapsis, deg 

Subscripts: 


1,2 

denotes specific points on the ellipse 

max 

maximum 

min 

minimum 


Symbols without arrows denote magnitudes. 

MISSION-ANALYSIS CAPABILITIES 

The computer program has been developed to fulfill a requirement for preliminary 
mission analysis and design. The program is intended to be a rapid engineering tool 
which may be easily modified to perform additional tasks as the need arises. In this 
light, the following paragraphs describe the basic assumptions and approximations, the 
method of calculation, and the program capabilities for each part of the mission-analysis 
problem treated by the program. 


Heliocentric Trajectory 

The Earth- Mars-Sun geometry used for calculating the heliocentric trajectory ele- 
ments is shown in figure 1 . Point masses and Keplerian mechanics are assumed through- 
out the analysis. The heliocentric orbits of Earth and Mars are represented by time- 
varying mean orbital elements. If the position vector to the Earth at a launch date and 
the position vector to Mars at an arrival date (fig. 1 ) are known, a number of methods can 
be used to generate a unique set of trajectory elements which connect these two points in 
the desired trip time. A true anomaly iteration method (ref. 1 ) is used here. Once the 
elements of the heliocentric transfer trajectory are known, many additional parameters of 
interest are computed. For example, C3 (twice the injection energy of the spacecraft 
relative to the Earth) and (the declination of the launch asymptote relative to the 

Earth’s equator) are calculated. These two parameters are of interest to the mission 
planner because they define launch- vehicle energy requirements per unit mass (C3) and 
whether the injection into the interplanetary trajectory violates range-safety requirements 
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Mars 


Sun 



Hel i ocent r I c 
transfer trajectory 


Earth 


Figure 1.- Earth to Mars geometry. 

(overflight restrictions on Dla)- Also, the hyperbolic excess velocity of the spacecraft 
relative to Mars Voo is computed here. Constraints on the maximum values of C 3 , 
Dla> applied by the program. The trajectory for a given launch and arrival 

date is rejected if it violates any one of the constraints, and a new launch-arrival date 
pair is tried. Thus, the mission planner is spared the necessity of sifting through a num- 
ber of impractical trajectories. The other parameters computed here are described in 
appendix C. This trajectory computation is very rapid and may be easily modified to 
compute additional quantities of interest to the mission planner. 


Elliptical Orbit at Mars 

The orbit trace and inertial landing-point geometry at Mars is shown in figure 2. 
The mission requirements of Sun lighting angle G, landing latitude X^, orbital inclina- 
tion i, and the argument of the landing point f^ are specified. The declination and right 
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Figure 2.- Landing -point geometry. Arrows indicate positive sense. 

ascension of the subsolar point, Xg and 0g, are calculated from the position of the Sun 
with respect to Mars. Since these quantities are known, the argument of periapsis co 
and the right ascension of the ascending node O can be found. (See fig. 2.) The 
resulting equations depend on the location of the landing point with respect to the ascendii^ 
or descending node of the orbit and with respect to the mornii^ or evening terminator 
(that is, lighting conditions). The various combinations of landing-point conditions are 
chosen on option from the main program. Finally, the apoapsis and periapsis radii, 
r^ and rp, are specified from experiment considerations. Thus, the orbital elements 
(ra, rp, i, J2, and co) of an ellipse which passes over the inertial landing point are 
known for the date of deorbit. 
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For photographic coverage, the spacecraft will be in orbit about Mars for a number 
of revolutions prior to landing. Because of the oblateness of Mars, and w will 
change as functions of time. Therefore, S2 and w are regressed an amount dependent 
upon the required stay time in orbit. (See ref. 1.) Thus, the elements of the initial 
ellipse on the date of the deboost maneuver are determined. 

The mission planner is interested in Sun and Earth occultations as seen by the 
spacecraft while in orbit about Mars. Therefore, such parameters as the first orbit on 
which occultation occurs, duration of occultation, and the time and true anomaly from 
perlapsis of entrance to and exit from occultation are computed for both the Sun and 
Earth. These parameters are necessary to define quantities such as battery require- 
ments (solar cells occulted from sunlight) and data-storage capability (direct telemetry 
occulted from tracking bases). The computed quantities are described in appendix C. 

It would be possible to compute occultation parameters for other celestial bodies (for 
example, Canopus) by a suitable modification to the program. 

Deboost Maneuver 

The deboost maneuver geometry is shown in figure 3. A minimum AV impulsive 
burn maneuver is computed. The maneuver is not constrained to be coplanar or to be a 
periapsis-to-periapsis transfer. The values of hyperbolic excess velocity Voo and a 
unit vector parallel to the approach asymptote and passing through the center of the 
planet S have been computed in the heliocentric trajectory part of the program. The 
quantities Voo and S define a family of approach h 3 rperbolas. The orbital elements 
of the required ellipse at Mars have been determined in the elliptical orbit computation 
section. The deboost maneuver is designed to specify the family of approach hyperbolas 
that results in the minimum AV requirement for de boost. 

The procedure is described with reference to figure 3. The approach hyperbola is 
rotated about S* and its periapsis altitude is adjusted until it intersects the specified 
elliptical orbit at a particular true anomaly f^. Since the radius of the hyperbola is 
constrained to be the same as the radius of the ellipse at that true anomaly, the orbital 
elements of the h 3 q>erbola are computed. The velocities on the ellipse Ve^l and h 3 q)er- 
bola Vj^^l at the intersection point are computed and their vector difference AVj 
represents an impulsive -burn transfer between the conics at their intersection. Next, 
another true anomaly f 2 is chosen, and the velocity difference AV£ is computed at this 
new intersection of the hsqjerbola and ellipse. This process is repeated at true anomaly 
intervals around the ellipse. The minimum AV is calculated by a parabolic interpola- 
tion through the three smallest computed values of AV. The associated hjqierbola is 
then defined to be the required conic. A maximum acceptable AV is defined by the 
user and any profiles which violate this constraint are rejected. 
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Figure 3-- Deboost velocity computation geometry. 


Additional parameters of interest to the mission planner are computed in this part 
of the program. For example, the areocentric components of AV, the plane change 
involved in the maneuver, and the radius, time and true anomaly of the deboost maneuver 
are computed. Also computed are the components B • T, B • R of a "miss parameter" 
B which is the perpendicular from the center of the planet to the approach asymptote. 
The components B • T and B • R lie in the plane formed by T (a unit vector perpen- 
dicular to S and parallel to the ecliptic plane) and R (r = x T, with R a unit vec- 
tor completing the triad). The RST areocentric coordinate system is a convenient 
targeting system for the mission planner. Other computed parameters associated with 
the deboost maneuver are described in appendix C. 
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Operational Modes 


There are several modes of operation and program options available to the mission 
planner. There are three output modes which control computational flow as shown in the 
flow chart in appendix A. A sample input and a sample output for each of the computa- 
tional modes are illustrated. An initial launch date of August 9, 1973, and an initial 
arrival date of March 16, 1974, are specified for each example. In each mode, the pro- 
gram automatically scans a grid of launch and arrival dates as determined by the user. 
Maximum values for the C3 , Dl^, Voo, and AV constraints are selected. A set of 
landing-point parameters (fj, X^, G, i, and stay time in orbit) are chosen which relate 
to the particular mission. Physical constants associated with the planet are specified. 

For each case, a set of program control integers is required. Each of the input and out- 
put parameters is described in appendixes B and C. 

An example of the input and output for the minimum output mode is presented. This 
operational mode performs only the calculations required to define the basic mission pro- 
file. Output is restricted to a single line to facilitate scanning a wide range of possible 
launch and arrival date combinations. Only the profiles which satisfy constraints on C3 , 
Dl^, Voo, are printed out. This option requires less than 1 second of computer 

time per launch- arrival date pair. 
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A sample input and output described in appendixes B and C for the minimum output 
mode follows; 


STNPT 


E 

= 

0.73E+02^ 0-86^-01* 0.9E + 01* O.Ot 

0.0* O.Ot 

XM 


0,74E + 02» 0,3E+01* 0,16E+02f 0.0 

♦ O.Ot 0.0* 

TLD 

= 

25. 


TAD 

= 

25* 


IJO 


5* 

ti 

0.428284E+05* 

J JO 


5* 

INC rr 

1 * 

C3MAX 

= 

0.21E+02* 

KEYl 

0* 

DLAMAX 

— 

0.4E+02. 

KEY3 

1* 

VHPMAX 


0.35E+01 . 

KEY4 = 

0* 

061 VMAX 

— 

0.2E+01* 

KFY5 

0* 

PFR 


0.12E + 02 * 

ISEARCH = 

0* 

ELP 

— 

0.3E+02. 

TOEST 

0.0* 

<;p6rL0N 

- 

0. 9086646901346 lF+02* 

iENO 


GEE 


C.65E+02 • 


XT 

- 

0.4F+02. 


DAYS 


0. lE+02* 


HA 


0.3267E+05 % 


HP 


0.14E+04, 


R S 

= 

0. 33934E+04* 


RPMIN 

r= 

0.47969F+04* 


T AOFDRB 

= 

0.0* 


FPA 

= 

-0.16E + 02 * 


RFNTRY 

= 

0.363724E+04* 


KJ? 

— 

0.197F-02* 
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LAUNCH 

ARRIVAL 

C3 

DLA OELTAV F DPST 

F drst drst 

TIM SO EO TSUNIN 

O'JRSOM 

TASTN TASPMT 

TFARIN 

-\(|Qr J5C 

TA-TA' T ACPI IT 

DATF 

DATE 


ELIPSF 

HYPERR FROM 

PFR 







73 

8 

9 

7A 

3 

16 

16.54 

33.66 

1.510 -63.7 -44.6 

-1630.69 

0 

1 

0.00 

0.00 

0.00 

0.00 

8.84 

'‘6. 11 

?4.5’ 

75.07 

73 

8 

9 

74 

3 

21 

17.04 

35.96 

1.689 -74.4 -51.1 

-2065.39 

0 

1 

0. 00 

0.00 

0.00 

0.00 

«.9? 

’7.21 

’4.71 

76.58 

73 

8 

9 

74 

3 

26 

17.71 

38.51 

1.877 -84.1 -55.5 

-2556.29 

0 

1 

0.00 

0.00 

0.00 

0.00 

9.0? 

7f\, 70 

’4.98 

78.05 


73 

8 

14 

74 

3 

16 

17.41 

28.85 

1.349 -60.8 -41.6 

-1529.65 

0 

1 

0.00 

0.00 

0.00 

0.00 

8.84 

26.11 

24.5’ 

75.07 

73 

8 

14 

74 

3 

21 

17.67 

30.67 

1.501 -71.5 -48.4 

-1939.69 

0 

1 

0.00 

0.00 

0.00 

0.00 

8.0? 

’7.21 

?4.’l 

76. 5S 

73 

8 

14 

74 

3 

26 

18.03 

32.68 

1.658 -81.3 -53.7 

-2404.05 

0 

1 

0.00 

0.00 

0.00 

0.00 

9.0? 

?B . ’o 

24.98 

78.05 

73 

8 

14 

74 

3 

31 

18.52 

34.90 

1.820 -90.2 -57.0 

-2925.13 

0 

1 

0.00 

0.00 

0.00 

0.00 

9. ^ ^ 

’9. ’5 

25. -*2 

79.49 

73 

8 

14 

74 

4 

5 

19.18 

37.35 

1.990 -98.0 -58.4 

-3502.91 

0 

1 

0.00 

0.00 

0.00 

0.00 

9.30 

30. ’9 

?5.7? 

BO. 88 


73 

8 

19 

74 

3 

16 

19.15 

24.76 

1.237 -59.4 -39.1 

-1481.56 

0 

1 

0.00 

0.00 

0.00 

0.00 

8.R4 

?5. 1 1 

?4.52 

75.07 

73 

8 

19 

74 

3 

21 

19.25 

26.22 

1.366 -69.8 -45.9 

-1870.06 

0 

1 

0.00 

0.00 

0.00 

0.00 

8.92 

? 7 . ? ] 

?4.71 

76.58 

73 

8 

19 

74 

3 

26 

19.41 

27.81 

1.502 -79.5 -51.4 

-2309.39 

0 

1 

0.00 

0.00 

0.00 

0.00 

Q.O? 

28. ?9 

’4.98 

78.05 

73 

8 

19 

74 

3 

31 

19.64 

29.56 

1.640 -88.2 -55.3 

-2801.39 

0 

1 

0.00 

0.00 

0.00 

0.00 

9.15 

’9. ’5 

?5.’? 

79.49 

73 

8 

19 

74 

4 

5 

19.97 

31.48 

1.781 -96.1 -57.5 

-3347.84 

0 

1 

0.00 

0.00 

0.00 

0.00 

9.30 

’0.39 

25.7? 

80. 88 



An example of the input and output for the extended printout mode is presented. 
This operational mode calculates numerous additional parameters (described in appen- 
dix C) associated with a particular launch and arrival date. In this mode, the "time- 
correction" option can be selected. This option allows the user to pick a particular lon- 
gitude of landing (that is, "tyir^ down" the inertial landing point to the rotating planet). 
The chosen longitude must rotate beneath the computed elliptic orbit with the given 
lighting conditions at a particular time. This constraint determines the landing time. 
With the time of landing known, event times of deorbit and deboost are computed. The 
deorbit -to-landing time increment is computed by use of an entry trajectory with no 
atmosphere to allow rapid computation. On option, a more accurate time increment can 
be put into the program. The deboost -to-deorbit time increment is computed aloi^ the 
specified elliptical orbit. One page of output and approximately 2 seconds of computer 
time per launch-arrival date pair are required in the extended output mode. Other 
quantities of interest to the mission planner could be computed and inserted into this 
output. 
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A sample input and output for the extended output mode (described in appendixes B 
and C) follow: 


SINPT 


E 

=5 

0* 73E + 02^ 0.8E+01 ♦ 

0.9F+01 t O.Ot 0*0, 0.0 . 

XM 

- 

0.74E + 02f 0.3E+01 , 

0*l6E+02t 0*0t O.Ot O.Ot 

ILO 

= 



lAD 

= 

It 


IJO 

= 

It 


jjn 

= 

It 


C3MAX 


0.21E+02t 


OLAMAX 

= 

0.4E+02t 

IJ = 0.428284E + 05 . 

VHPMAX 


0.35E+01 . 

INC = It 

DELVMAX 

= 

0.2E^-01* 

KEYl = It 

PER 

= 

0.12E+02 . 

KFY3 = 1* 

ELP 

= 

0.2E+02t 

KEY4 = 1* 

SPECLON 

= 

0.32E+03t 

KFY5 = Ot 

GEE 


0 *6 F +02 1 

ISEARCH = Ot 

XI 


0.3E+0?t 

TDFST = 0.0* 

DAYS 

= 

0.9E+01 t 

^END 

MA 

= 

0.3267F+05t 


HP 


0* 14E+04t 


RS 


0*33934E+04t 


RPMIH 

= 

0.47969E+04t 


TADEORB 


0.21475E+03t 


FPA 

= 

-0.16E+02t 


RENTRY 


0.363724E+04t 


XJ2 

= 

0.197E-02t 
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LAUNCH DATE AR«M \/AL OATP r,cnqn'^T TIMF nrromT rpir lANOT^r, TTm- 

CAIENOAR 8 q 73 0 0 0 3 1 6 74 0 0 0 3 16 74 19 49 4 3 ?9 74 1 7 4'> ^ 74 ?0 18 4B 

JULIAN 2441903.50 2442132.50 2442123.3’ 744?1’3.?1 74421’?. ’5 

INTERPLANETARY FLIGHT PARAMETERS 

OLA = 3.36573132E + Q1 RAL = 1 . 3656 72 55F + 01 C’ = 1 . 65407 750: +01 Tcrp TTmf = 7.190000 00*^ + 0’ 

ARFO. OFC. S-VECTHR = -1 .08404062F+01 AREn. u.A, S-ViCTOR = 1 . 7 1 786 355F +0 2 HYPfR. vfl , ^ 7 . 5’ 6P ’ 0’ T*" + 00 

GFO. nrC. S-VFCTOR = - 1.74.848575E +01 0^0. R.A. S-VRCTHP = 4 . 1 7 90 7 72*’ E +0 1 roM T PA t j mm Ot<;t, = ’ . ^006 PT*" +0B 

7AP = 1 .03687925F+0? ET<; = 1 . 73 5583 82F +02 2 AE = 1.32057526^ + 02 " = 2 . 00 ” 1 ’ 46*' +0’ 

7AC = 5.4663C549F+01 ETC ^ 2.694755439+02 

PROBE PERIHELION = 1 . 51 3 1 81 94 E+OB PROBE APHELION = 2.457373659+08 ORO^r inOLTmattom = ’.05159177^+00 

LAUNCH TRUE ANOM . s 8 . 84 53 1 1 BOE+ 00 ARRIVAL TRUE ANOm. = 1 . 5 8446 95 1 *■ +o 2 H«=Ln, ANGL^ tRAVCi = i .49601 6 39F + 0 7 

R-VECTCR MAGNITUDE = 9 .8 3? 10 86 3^ +Q 3 8 OUT T = 9. 5 853? 051 F + 03 « 0'’T R = -’.18)161709+0’ 

fclFMENTS AND DEROOST PARAMETERS - HYREPROLA 

A = -6.65501483E+03 E = 1 . 78 40 14C7E+ 00 I = 3. 16 568 2 36E +01 CAR.OMf OA s -7 .6 ’ 0746 ’1 F +01 = 1 . 450 91 7’ ?? +02 

DRST TRUE ANOM. s -3. 73 71 3246E + 01 OPST TIMF = - 7 . 8 3871 ? 3 1 E+02 

PFR. RADIUS =: 5.21762525F+03 PER. = 1 . 74781230E + 01 P^R. R.A. = 7. 79875332F+01 

V AT ORST = 4.54892991E + 00 VX AT pnST = -3 . 3 1 1 9” 1 9E + 00 VY AT ORST = 7. 6808665 59 + 00 V7 A’’ 0 ‘’*:t - - 1 . 592 7? 860'’ +00 

ELEMENTS AND OEBOOST PARAMETERS - ELLIPSE 

A = 2.04284000E + 04 E = 7 .65 35607 3E-01 ! = 3.00000000E +01 CAP.OMi-OA = 3. 05 3851 71 F + 07 OMcr.^ 

DRST TRUE ANOM. = -5. 7743 5330 E +0 1 ORST TIM9 = - 1 .4 249 8943E + 03 DORR, T. A. = 2 . 1 4 7500 OOF +0 ’ vnncQ 

PER. RADIUS = 4. 79340000E+03 PER. O^C. = 1 . 589982 17F+0 1 '’Fp. r.a. = 9.582206859+01 orotno 

V AT DRST = 3.48722587E+00 VX AT DRST = - 2 . 83 902 1 72E + 00 VY AT DRST = 1 . 899 743’ 4F +00 V7 AT ORCt 

DFBOnST MANEUVER PARAMETERS 

DELTA V = 1. 27618493F + 00 EXCESS DELTA V = 7. 2 38 1 5067F-0 1 RADIUS = 6 . 007943R9F + 03 opo = I , ?0000000*- + 01 

VX = 4.72911478F-01 VY = -7 . 8 1 1 2321 OE-01 V7 = 8 . 91 543065E-0 1 pLAN^ = 1 , 11 8 74’ 889 + 0 1 

nCCULTATICN PARAMETERS 

1ST ORBIT, SUN= 0 TIME, SUN = 0. DiJRATION, SUN = 0. 

TRUE ANOM., SUN IN = 0. TRUE ANQM., SUN OUT = 0. 

1ST ORBIT, EARTH = 1 TIME, EARTH = 1 . 25 7587 2 lE + 01 DURATION, EARTH = ? .44956743F +0 1 

TRUE ANOM., EARTH IN = 3 -3997903 3E+0 1 TRUE ANOM., EARTH OUT - 7 . 7 764611 5E +01 

LANDING POINT PARAMETERS 

LCNGITUDE = 3 .2 OOOOOOOE +0? LATITUDE = ? .000 OOOOOE +0 1 SUN LIGHTING ANGLE 6.000000009 + 01 nAY*^ = 9 


= 1.46776084=+02 
= 2.?0634Pl’9-oi 
= 1 .02581 6 ’49 +00 
= -7.011955689-01 



The third output option is the "parametric-analysis" printout mode. An example 
of the input and output of this mode follows. This mode computes the heliocentric trans- 
fer trajectory for a specified launch and arrival date. Then the input parameters, 

Xi, G, i, and stay time in orbit, are varied through specified raises to determine their 
effect on the AV required for the deboost maneuver. One line of output and approxi- 
mately 1 second of computer time are required for each combination of input parameters. 


SINPT 



E 

= 

0.73E+02f 0 

XM 


0.74E+C2* 0 

TLO 



lAD 

= 

1 » 

un 

= 

1. 

J JO 

= 

It 

C3MAX 

= 

0.21E^02 t 

DLAMAX 

= 

0.4F+02t 

VHPMAX 

= 

0. 35E^-01 * 

OELVMAX 


C.2E+01 . 

PFR 

= 

0.12E+02t 

El P 


0.2E+02t 

SPECl. ON 

= 

0.32F+03. 

OFE 

= 

0.6F+02t 

XI 


0* 3F + 02t 

DAYS 


0.9E+01* 

HA 


0.3267E+05, 

HP 


0. 14E+04t 

RS 


0. 33934E+04. 

RPMIN 

= 

0.47969E+04* 

TAOEORB 

= 

0.21475F+03, 

FPA 

= 

-0.16E+02* 

RENTRY 

= 

0.363724E+04 

XJ2 

- 

0.197E-02t 
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tj 


0*4282 84E+-05 

INC 


1 * 

KFYl 


?♦ 

KFY3 


1. 

KEY4 

= 

1 * 

KFY5 

= 

0^ 

T SEARCH 


1 ♦ 

THEST 

.= 

o 

• 

o 

SFNO 



iPARAM 



PEP 1 

= 

0. lE+02* 

PER? 


0*12E+02* 

KPFR 

= 

2* 

FLPl 


0*25E+02, 

ELP? 


0.3E+02, 

KELP 


5, 

GFFl 


0.6F+02» 

GEE2 

= 

0*65F+02* 

KGEF 


5t 

XTl 

= 

0.356+02, 

XI? 

= 

0. 46+02, 

KXT 

= 

5, 

DAYl 

= 

0.5E+01, 

DAY? 

= 

O.lF+02, 

KHAY 


5, 

SFNO 
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CALENDAR 
MILT AN 


LAUNCH DATE 
8 9 73 0 0 0 

24A1903.50 


ARRIVAL DATE 
3 16 74 0 0 0 

2442122.50 


PER 

S 

1.0000E4-01 

FLP 

= 

2. 5000E<-01 

GEE 

= 

6.0000E+01 

PER 

= 

l.OCOOE+01 

ELP 


2. 5000E+01 

GEE 


6.0000E+01 

PER 

= 

l.OOOOE+01 

ELP 


2. 5000F4-01 

GEE 

= 

6.0000E+01 

PER 

« 

l.OOOOE+01 

ELP 

= 

2.5000E+01 

GEE 

= 

6.0000E+01 

PER 

X 

l.OOOOE+01 

ELP 

= 

2.5000F+01 

GEF 


6.5000E+01 

PER 

s 

l.0000E<'01 

ELP 

3 

2,50036+01 

GEE 

= 

6.5000E+01 

PER 

X 

1.0000E4-01 

ELP 

= 

2,5000F+01 

GEE 

= 

6.5000E+01 

PER 

3 

l.OOOOE+01 

FLP 


2.50006+01 

GEE 

= 

6.5300E+01 

PER 

= 

l.OOOOE+01 

ELP 

= 

3.00036+01 

GEE 

= 

6.0000E+01 

PER 

3 

l.OOOOE+01 

FLP 

= 

3.0000E+01 

GEE 

= 

6.0000E+01 

PER 

3 

l.OOOOE+Ol 

FLP 


3.0003E+01 

GEF 

= 

6.0000E+01 

PER 

3 

l.OOOOE+01 

ELP 

= 

3.0000E+01 

GEE 

= 

6.0000E+01 

PER 

X 

l.OOOOE+01 

ELP 

= 

3.0000E+01 

GEF 

= 

6.5000E+01 

PER 

3 

l.OOOOE+01 

ELP 


3.00006+01 

GEF 

= 

6.5000E+01 

PER 

3 

l.OOOOE+01 

FLP 

= 

3.0000F+01 

GEE 


6.5000E+01 

PER 

3 

1.0000E4-01 

ELP 

= 

3.0000E+01 

GFP 

= 

6. 5000E+01 

PER 

3 

1.2000E+01 

ELP 

= 

2.5000E+01 

GEE 

= 

6. OOOOE+01 

PER 

3 

1,2000F+01 

FLP 


2.5000E+01 

GEE 

= 

6. OOOOE+01 

PER 

3 

l-POOOE+Ol 

ELP 


2.5000E+01 

GEE 


6.0000E+01 


-4 


XI 


3.5000F+01 

DAYS 

= 

5.0000E+00 

DELTA 

V 


1.2627E+00 

XI 

= 

3.5000E +01 

DAYS 


l.OOOOE+01 

DELTA 

V 

= 

1.3353E+00 

X I 

= 

4.00006+01 

DAYS 

= 

5.0000F+00 

DELTA 

V 

= 

1.0735E+00 

XI 

- 

4.0000E+01 

DAYS 

= 

l.OOOOE+Ol 

DELTA 

V 

= 

1.1385E+00 

XI 


3.50006+01 

DA Y^ 

= 

5.0000E+00 

DELTA 

V 


l.3832E+0n 

XI 


3.5000F+01 

DAYS 

= 

l.OOOOE+Ol 

delta 

V 


1 .4547B+00 

XI 

= 

4.0000F+01 

DAYS 

= 

5.0000E+00 

delta 

V 


1.1932F+00 

XI 

= 

4.0000F+01 

DAYS 

= 

l.OOOOE+Ol 

DELTA 

V 

= 

1.2675F+00 

XI 

= 

3.5000«= + 01 

DAYS 


5.0000E+00 

DELTA 

V 


1.6100E+00 

XI 

= 

3.5000E+01 

DAYS 


l.OOOOE+Ol 

DELTA 

V 

= 

1.6812E+00 

XI 

= 

4.0000E+01 

DAYS 

- 

5.0000E+00 

DELTA 

V 

3 

1.2901E+00 

XI 

= 

4.0000E+01 

DAYS 


l.OOOOE+Ol 

delta 

V 

= 

1.3765E+00 

XI 

= 

3.50006+01 

DAYS 


5.0000E+00 

DELTA 

V 


1.7168F+00 

XI 

= 

3.5000F+01 

DAYS 


1 .OOOOE+Ol 

DELTA 

V 

= 

1.78156+00 

XI 

= 

4.0000F+01 

DAYS 

= 

5.0000E+00 

delta 

V 

= 

1.4342F+00 

XI 

= 

4.0000E+01 

PAYS 

= 

l.OOOOF+Ol 

DELTA 

V 


1.5179E+00 

XI 


3.5000E+01 

DAYS 


5.0000E+00 

DELTA 

V 

3 

1.2611E+00 

XI 

= 

3.5000E+01 

DAYS 

= 

l.OOOOF+Ol 

delta 

V 

= 

1.3326F+00 


XI 


4.0000E+01 


DAYS 


5.0000F+00 DEL'TA V 


1.0812E+00 



There are several program options available within each of the output modes. The 
program user must select values for each of the control integers described in appen- 
dix B. Either a posigrade or retrograde hyperbola may be chosen for the approach to 
Mars. One of four combinations of ascending or descending node and morning or evening 
lighting are specified for the landing point on the surface of Mars. The user may also 
put in broken-plane trajectory parameters (that is, a combination of two intersecting 
conic trajectories with different orbital elements). If this input option is selected, no 
heliocentric trajectory computation is made within the program. (See the flow chart in 
appendix A.) 

A brief description of the primary subroutines is given in appendix A. The other 
general purpose subroutines are given in the FORTRAN listing in appendix D, These 
subroutines are generally self-explanatory and may be used in many ts^jes of programs. 
The mathematics associated with several of the general purpose subroutines is discussed 
in reference 2 . 


Application to Other Planets 

The modular construction of the mission-analysis program permits a relatively 
easy extension to the evaluation of a planet -to-planet mission. A planetary ephemeris 
subroutine must be substituted for the launch and arrival planets, and a coordinate trans- 
formation subroutine is required to rotate a vector between the mean planet equator- 
equinox coordinate systems. Changes in the planetary constants and minor FORTRAN 
modifications must also be made. The straightforward computation and flexibility of the 
mission-analysis program should permit the addition of any calculations required by the 
mission planner. 


RESULTS AND DISCUSSION OF SAMPLE CASE 

Many phases of preliminary mission analysis can be studied by use of the data com- 
puted by the program. Since this paper is intended as an explanation of the capabilities 
of the mission-analysis program, only two examples of data analysis are presented here. 
First, consider the problem of choosing feasible ranges of launch and arrival dates. 

This analysis is performed by the minimal output mode of the program. The choice of 
a launch and arrival date pair is constrained immediately by the available launch vehicle 
energy per unit mass C 3 , range safety considerations D^^, and the deboost velocity- 
change requirement AV. Program output is used to plot constant contours of these three 
quantities as functions of launch and arrival dates. (See fig. 4.) By establishing an upper 
value of each quantity, many combinations of launch- arrival dates which exceed one of the 
constraints can be eliminated immediately from further consideration. The shaded region 
in figure 4 indicates the launch and arrival dates which simultaneously satisfy the 
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Orbital inclination, deg 




constraints of C 3 = 20 kmVsec^, = 35°, and AV = 1.3 km/sec. It should be 

noted that the AV contours are dependent upon a particular set of landir^-point condi- 
tions which determine the elliptical orbit at Mars. Any variation in the elliptical orbit 
will shift the AV contours. 

Next, it is logical to choose a launch-arrival opportunity from the shaded region 
of figure 4, and to determine the effect of the landing-point conditions on AV. This 
problem is easily handled by the parametric analysis output mode. Figure 5 shows the 
variation in inclination as a function of deboost AV for a launch date of August 20, 1973, 
and an arrival date of March 24, 1974. The plots show the variation in AV as depen- 
dent upon the Sun lighting angle, ii, landing latitude, and stay time in orbit. Since 
these parameters are not truly independent, the two-dimensional plot will not tell the 
complete story. However, this type of analysis does indicate trends for additional study. 
The data presented in figures 4 and 5 can be generated in two runs of the program. 

CONCLUDING REMARKS 

Rapid computing time, modular construction, and the combination of many phases 
of mission design into one package are assets of the mission- analysis program. The 
conic trajectory computations decrease the accuracy of the program as compared with 
an integrated trajectory, but this inaccuracy is felt to be of minor concern for prelimi- 
nary mission planning purposes. The ease with which modifications can be made, and 
computational flexibility make the program a useful engineering tool. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 7, 1970. 
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APPENDIX A 


MAIN PROGRAM 
Flow Diagram 



This flow diagram is a simple description of the computa- 
tional flow of the main program. The three sets of input are 
described in appendix B. Each of the subroutines shown here 
is described briefly after the flow diagram. These subroutines, 
in turn, call the other subroutines given in appendix D. The 
numbered decisions are as follows: 

Is the program in the parametric analysis output mode ? 
(2) Is broken-plane input data required? 

Is the program in the extended printout mode? 

Are the C3, Dl^, or Voo constraints exceeded? 

(5) Is the AV constraint exceeded? 

Has the "time correction" option been selected? 

Is the required range of launch -arrival dates completed? 
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APPENDIX A - Concluded 


Description of Primary Subroutines 

The following is a brief description of the primary subroutines shown in the flow 
diagram. 

PLUG and PLUG 1 — PLUG is called in the minimum output mode; it computes C 3 , 
Dla> Voo, and the declination and right ascension of S . PLUG 1 is called in the 
extended output mode; it computes many additional trajectory parameters as described 
in appendix C. These subroutines compute the elements of a heliocentric trajectory 
between Earth and Mars for a given launch and arrival date. A unique conic trajectory 
is determined from the two radius vectors and the trip time. A true anomaly iteration 
method is used to establish the conic trajectory. 

COMPEL — If the landing-point parameters of i, f^, G, X^, stay time in orbit, and the 
date of landing are known, this subroutine computes the elements of an ellipse at Mars 
which passes over the landing point on the date of deorbit. 

FUDGE — With the oblateness coefficient for Mars and the stay time in orbit known, this 
subroutine modifies and u> as computed in COMPEL to account for oblateness 
effects accumulated during the specified stay time in orbit. The resultant elements 
determine the initial orbit at Mars. 

DEBOOST AND DBST 1 - DEBOOST is called in the minimum output mode; it computes 
the minimum AV and the elements of the hyperbola associated with it, and the impact 
plane parameters. DBST 1 is called in the extended output mode; it computes the addi- 
tional deboost parameters described in appendix C. These subroutines compute the 
minimum impulsive velocity change required to transfer from the approach hyperbola to 
the initial orbit at Mars. 

SEESE — This subroutine determines whether occultations of the Sun or Earth as seen 
by the spacecraft occur during the specified stay time in orbit. If so, it computes the 
occupation parameters described in appendix C. 

EVENTS — This subroutine is called in the "time correction" option. The user specifies 
a longitude of landing which determines the landing time. Then, event times of deorbit 
and deboost are computed. The user may specify the deorbit -to-landing time increment. 
If this time increment is not specified, an impulsive deorbit maneuver is performed and 
the deorbit quantities described in appendix C are computed. 
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APPENDIX B 


DESCRIPTION OF PROGRAM INPUT 


The NAMELIST feature of the Control Data 6600 computer is used to facilitate data 
input. Three namelists are used. INPT is the standard set of data and control variables 
and is always required. PARAM is a set of quantities used to vary the landing point data 
parametrically and is put in only in the parametric-analysis output mode. BPDATA is 
a set of broken-plane trajectory parameters. If this namelist is put in, no trajectory 
calculation is made within the program. 

The following quantities are put in by the INPT namelist. 


Program Program Mathematical 
symbol dimension symbol 


Units 


Definition 


E 

6 


yr, mo, day, 
hr, min, sec 

Calendar launch date (for example, 73, 
10, 23, 0, 0, 0) 

XM 

6 


yr, mo, day, 
hr, min, sec 

Calendar arrival date 

ILD 

1 


days 

Scanning period size for launch 

lAD 

1 


days 

Scanning period size for arrival 

UD 

1 


days 

Launch-date increment 

JJD 

1 


days 

Arrival-date increment 

C3MAX 

1 

1 

^ 3 , max 

kmVsec^ 

Constraint on vis-viva injection 
energy C 3 

DLAMAX 

1 

^LA,max 

deg 

Constraint on declination of the launch 
asymptote Dj__^ 

VHPMAX 

1 

Voo,max 

km/sec 

Constraint on h 3 q)erbolic excess 
velocity Voo 

DELVMAX 

1 

^^max 

km/sec 

Constraint on deboost AV 

PER 

1 

h 

deg 

Periapsis to landing point angle 

ELP 

1 

h 

deg 

Landing-point latitude 

SPECLON 

1 

<^p 

deg 

Specified landing-point longitude 
(required only in time-correction 
mode) 

GEE 

1 

G 

deg 

Sun lighting angle 
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APPENDIX B - Continued 


Program 

S3rmbol 

Program 

dimension 

Mathematical 

symbol 

Units 

Definition 

XI 

1 

i 

deg 

Orbit inclination at Mars 

DAYS 

1 


days 

Stay time in orbit 

HA 

1 

^a 

km 

Height of apoapsis above Mars surface 

HP 

1 

hp 

km 

Height of periapsis above Mars surface 

RS 

1 

V 

km 

Radius of Mars 

RPMIN 

1 

^p,min 

km 

Minimum periapsis radius of approach 
h3rperbola 

TADEORB 

1 

f 

deorbit 

deg 

True anomaly of deorbit (required only in 
time-correction mode) 

FPA 

1 

a 

deg 

Flight -path angle at entry 

RENTRY 

1 

^entry 

km 

Radius at entry (FPA and RENTRY are 
required only in time -correction mode 
and when no estimated time from 
deorbit to landing has been specified) 

XJ2 

1 

J2 

none 

Oblateness coefficient for Mars 

U 

1 


km^/sec^ Gravitational constant for Mars 

INC 

1 


none 

INC is 1 for posigrade h 3 q>erbola and 2 
for retrograde h 3 q)erbola 

KEYl 

1 


none 

Output control integer 

KEYl is 0 for minimum-output mode, 

1 for extended-output mode, and 

2 for parametric-analysis output mode 

KEYS 

1 


none 

Landing-point control integer 
KEYS is 1 for descending node, 
p.m. lightii^; 2 for ascending node, 
p.m. lighting; 3 for descending node, 
a.m. lighting; 4 for ascending node, 
a.m. lighting 

KEY4 

1 


none 

Time- correction mode control integer 
KEY4 is 1 for time correction, and 
0 for standard run 
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APPENDIX B - Concluded 


Program Program Mathematical units 
symbol dimension symbol 


Definition 


KEYS 1 none Broken-plane input control integer 

KEYS is 1 for broken-plane input (BPDATA 
namelist must be added) and 0 for standard 
input 

ISEARCH 1 none Parametric-analysis control integer 

ISEARCH is 1 for parametric analysis (PARAM 
namelist must be added) and 0 for standard run 

TDEST 1 days Estimated time from deorbit to landing 

TDEST is 0 yields computed time from deorbit 
to landing 

The following quantities are input by the PARAM namelist when required. 


Program symbol 

Units 

Definition 

PERI 

deg 

Initial value of f^ 

PER2 

deg 

Final value of f^ 

KPER 

deg 

Incremental value of ^ 

ELP1,ELP2,KELP 

deg 

Initial, final, and incremental values of 

GEE1,GEE2,KGEE 

deg 

Initial, final, and incremental values of G 

XIipa2,KXI 

deg 

Initial, final, and incremental values of i 

DAY1,DAY2,KDAY 

days 

Initial, final, and incremental values of stay time 



in orbit 

The following quantities are input 

by the BPDATA namelist when required. All are 

associated with a particular broken-plane trajectory. 

Program symbol 

Units 

Definition 

C3 

km2/sec2 

C3 

DLA 

deg 

DlA 

VHP 

km/sec 

Voo 

DPA 

deg 

Declination of the approach asymptote 

RAP 

deg 

Right ascension of the approach asymptote 

XM 

yr, mo, day, hr, min 

, sec Arrival date for broken-plane trajectory 

E 

yr, mo, day, hr, min 

, sec Launch date for broken-plane trajectory 
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APPENDIX C 


DESCRIPTION OF PROGRAM OUTPUT 


Three output options are available in the program. The first output option is a 
minimum print mode. Several important parameters associated with a particular launch- 
arrival date pair are printed on a single line. This mode facilitates scanning a wide 
range of launch-arrival date combinations to select suitable mission profiles. Only pro- 
files which satisfy the C 3 , Dl^, Voo, and AV deboost constraints are printed out. 

The program automatically scans a grid of launch and arrival dates as determined by the 
first six quantities in the INPT namelist. The following quantities are output in the mini- 
mum print mode. 


Output 

Units 

LAUNCH DATE 

yr^ mo, day 

ARRIVAL DATE 

yr, mo, day 

C3 

kmVsec^ 

DLA 

deg 

DELTAV 

km/sec 

F DBST ELLIPSE 

deg 

F DBST HYPERB 

deg 

DBST TIM FROM PER 

sec 

SO 

none 

EO 

none 


TSUNIN 

min 

DURSUN 

min 

TASm 

deg 

TASOUT 

deg 

TEARIN,DUREAR, 

min 

TAEIN,T ABOUT 

deg 


Definition 

Launch date at Earth 
Arrival date at Mars 

Vis-viva injection energy for Earth-Mars trajectory 

Declination of launch asymptote 

Deboost velocity change requirement 

True anomaly of deboost on the elliptical orbit 
at Mars 

True anomaly of deboost on the approach hjqjerbola 

Time of deboost from the periapsis of the ellipse 

The first orbit on which occultations of the Sun take 
place; SO = 0 indicates no occultations of the 
Sun during the specified stay time in orbit 

The first orbit of Earth occultations; EO = 0 

indicates no occultations during the specified stay 
time in orbit 

Time from elliptical periapsis of entrance to Sun 
occultation 

Duration of Sun occultation 
True anomaly of entrance to Sun occultation 
True anomaly of exit from Sun occultation 
Parameters associated with Earth occultations 
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APPENDIX C — Continued 


The second output option is an extended printout mode. This option performs the 
same taslra as the minimum print mode, but with many additional parameters computed. 
The program must be in this mode in order to select the time-correction option. This 
operational mode is useful for examinii^ a candidate mission profile in detail. The fol- 
lowing quantities are output in the extended printout mode. 


Output 

Units 

Definition 

LAUNCH DATE 


Calendar (mo, day, yr, hr, min, sec) and Julian 
(days) launch date from Earth 

ARRIVAL DATE 


Calendar and Julian arrival date at Mars 

DE BOOST TIME 


Calendar and Julian deboost date; output only in 
time-correction option 

DEORBIT TIME 


Calendar and Julian deorbit date; output only in 
time-correction option 

LANDING TIME 


Calendar and Julian landing date; output only in 
time -correction option 

DLA 

deg 

Geocentric declination of the launch asymptote 

RAL 

deg 

Geocentric right ascension of the launch asymptote 

C3 

km2/sec2 

Vis-viva injection energy 

TRIP TIME 

days 

Trip time 

AREO .DE C .S- VE CTOR 

deg 

— ► 

Areocentric declination of S 

AREO.R.A.S-VECTOR 

deg 

Areocentric right ascension of S 

HYPER.EXCESS VEL. 

km/sec 

Hyperbolic excess velocity 

GEO.DEC.S-VECTOR 

deg 

Geocentric declination of 

GEO. R.A.S- VECTOR 

deg 

Geocentric right ascension of S 

COMMUNICATION DIST. 

km 

Line-of-sight distance from Mars center to Earth 
center at arrival date 

ZAP 

deg 

Angle between S and Mars-to-Sun vector 

ETS 

deg 

Angle measured clockwise from T-axis to negative 
of projection of Mars-to-Sun vector on the RT 
plane (measured in areocentric, equatorial, 
arrival date coordinates) 

ZAE 

deg 

Same as ZAP with Mars-to-Earth vector 
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APPENDIX C - Continued 


Output 

Units 

ETE 

deg 

ZAC 

deg 

ETC 

deg 

PROBE PERIHELION 

km 

PROBE APHELION 

km 

PROBE INCLINATION 

deg 

LAUNCH TRUE ANOM. 

deg 

ARRIVAL TRUE ANOM. 

deg 

HELIO. angle travel 

deg 

B-VECTOR MAGNITUDE 

km 


B DOT T 

km 

B DOT R 

km 

The following parameters are output 
tical orbit about Mars. 

Output 

Units 

A 

km 

E 

none 

I 

deg 

CAP .OMEGA 

deg 

OMEGA 

deg 

DBST TRUE ANOM. 

deg 

DBST TIME 

sec 


Definition 

Same as ETS with Mars-to-Earth vector 

Same as ZAP with Mars-to-Canopus vector 

Same as ETS with Mars-to-Canopus vector 

Periapsis of heliocentric transfer trajectory 

Apoapsis of heliocentric transfer trajectory 

Inclination to the ecliptic of heliocentric 
transfer trajectory 

True anomaly of launch point on transfer 
trajectory 

True anomaly of arrival point on transfer 
trajectory 

Heliocentric angle between launch and 
arrival points 

Magnitude of B ("miss distance" from cen- 
ter of planet perpendicular to the approach 
asymptote) 

Component of B along the T-axis; areo- 
centric ecliptic of date coordinate system 

Component of B along the R-axis; areo- 
centric ecliptic of date coordinate system 

for both the approach h3rperbola and the ellip- 
Definition 

Semimajor axis of conic 

Eccentricity of conic 

Inclination to Martian equator of conic 

Right ascension of the ascending node of conic 

Argument of periapsis of conic 

True anomaly of deboost point on conic 

Deboost time from periapsis on conic 
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APPENDIX C - Continued 


Output 

Units 

Definition 

PER. RADIUS 

km 

Periapsis radius of conic 

PER. DEC. 

deg 

Areocentric declination of periapsis of conic 

PER. R.A. 

deg 

Areocentric right ascension of periapsis of 
conic 

V AT DBST 

km/sec 

Magnitude of velocity on conic at deboost 

VX,VY,VZ AT DBST 

km/sec 

Areocentric components of velocity on conic 
at deboost 

In addition to these parameters, 

the following quantities are output for the ellipse. 

Output 

Units 

Definition 

DORB. T.A. 

deg 

True anomaly of deorbit point on ellipse; 
output only in time-correction mode 

VDORB 

km/sec 

Velocity change required for deorbit; output 
only in time-correction mode 

PERIOD 

days 

Period of the ellipse 

The following parameters are output under the headings of "DEBOOST MANEUVER," 
"OCCULT ATION," and "LANDING POINT." 

Output 

Units 

Definition 

DELTA V 

km/sec 

Magnitude of velocity increment required for 
deboost 

EXCESS DELTA V 

km/sec 

^Vmax - ‘^V^eboost 

RADIUS 

km 

Radius on conics at deboost point 

PER 

deg 

True anomaly of landing point beneath the 
ellipse — (positive if landing before 
periapsis, negative if landing after 
periapsis) 

VX,VY,VZ 

km/sec 

Areocentric components of DELTAV 

PLANE CHANGE 

deg 

Angle between the planes of the approach 
hyperbola and the elliptical orbit 

1ST ORBIT, SUN 

none 

The first orbit on which occupations of the 
Sun occur 
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Output 

Units 

Definition 

TIME, SUN 

min 

Time of entrance to Sun occultation from 
periapsis 

DURATION, SUN 

min 

Duration of Sun occultation 

TRUE ANOM., SUN IN 

deg 

True anomaly of entrance to Sun occultation 

TRUE ANOM., SUN OUT 

deg 

True anomaly of exit from Sun occultation 

1ST ORBIT, EARTH 

none 

The first orbit on which occupations of the 
Earth occur 

TIME, EARTH 

min 

Parameters associated with Earth 

DURATION, EARTH 

min 

occupations 

TRUE ANOM., EARTH IN 

deg 


TRUE ANOM., EARTH OUT 

deg 


LONGITUDE 

deg 

Areocentric right ascension of the landing 
point (in time-correction mode, the speci- 
fied longitude is substituted) 

LATITUDE 

deg 

Specified latitude of the landing point 

SUN LIGHTING ANGLE 

deg 

Specified lighting angle at the landing point 

DAYS 

days 

Stay time in orbit prior to deorbit 


The third output option is a parametric-analysis printout mode. This mode lists 
the launch and arrival date for a particular trajectory. Then, on a single line, PER, 
ELP, GEE, XI, DAYS, and DELTA V are listed. Each landing-point parameter can be 
varied in turn to determine its effect of AV required for the deboost maneuver. 


oooooonnoooooooooooooooonooonoono 


APPENDIX D 


PROGRAM LISTING 


The following is a FORTRAN listing of the mission-analysis program and asso- 
ciated subroutines: 


PROGRAM MISHAP ( IN PUTt OUTPUT t TAPE5= I NPUT f TAPE6=0UTPUT ) 

DIMENSION E(6)>XM(6), El(6)tXMI(6I 

NAMELIST/INPT/E,XM, ILO. I AD , I JD , J JD i C3M AX , DL AMAX , VHPMA X , DEL VMAX , 
IPERfELP.SPECLONtGEE f X 1 1 DAYS . HA i HP , RS i RPMI N , TAOEOR B, FP A , RENTRY i X J2 t 
2Ut INCfKEYl ,KEY3,KEY4,KEY5, ISEARCHiTDEST 

4 /PARAM/P=Rl,PER2tKPcRt ELPl»ELP2tKELP»GEEltGEE2,KGEEf XI 1,X 

512 »KXI ,DAY1 tDAYZtXOAY 
2/BPDATA/C3,DLA,VHPt DPA,RAPtXM,E 
DIMENSION DRSTI6) ,DE0R(6) .XLAN0(6) 

1 READ(5,INPT> 

WRITEI6,INPT) 

IFdSEARCH.NE.DGO TO 11 
READIS.PARAM) 

WRITE(6,PARAM) 

11 CONTINUE 


INPUT FOR MAIN 

E(1-6)=FIRST CALENDAR CATE OF LAUNCH PERIOD. 

E( 1»=CALE^CAR YEAR12 DIGITS!. E(2)=CAL. MONTH, E(3)=CAL DAY, E(4«=H0URS 
E(5)=M1NUTES, E 16 )= SECONDS. 

XM( 1-6»=FIRST CALENDAR DATE OF ARRIVAL PERIOD. 

XMa,2,3,4,5,6)-SAME AS FOR LAUNCH. 

ILD=LAUNCH PERIOD SIZE, IAD=ARRIVAL PERIOD SIZE. 

IJD=LAUNCH DATE INCREMENT, JJD=ARRIVAL DATE INCREMENT, 

PER=IMPACT ANGLE, E LP = OECL INAT I ON OF IMPACT, SPECL ON = SP£C I F I ED 
LONGITUDE OF IMPACT, GEE=SUN LIGHTING ANGLE, XI = I NCL I NAT ION OF 
ELLIPTICAL ORBIT, HA=HEIGHT OF APOAPSIS, HP=HEIGHT OF PERIAPSIS, 
DAYS=STAY TIMS PRIOR TO DEORBIT. 

RS=RADIUS OF MARS, RPMIN=MINIMUM RADIUS OF PERIAPSIS OF HYPERBOLA, 
TADEORB=TRUE ANOMALY OF DEORBIT, FPA=FLIGHT PATH ANGLE AT ENTRY, 
RENTRY=RADIUS AT ENTRY, 

XJ2=06LATENESS COEFFICIENT, U*MU FOR MARS 

KEY1»0 FOR MINIMAL OUTPUT, =1 FOR EXTENDED OUTPUT, =2 FOR PARAMETRIC 
ANALYSIS OUTPUT. 

KEY3= 1, DESCENDING NOOE(PM) - 2, ASCENDING NQDEIPMI - 
3, DESCENDING NODEIAMI - 4, ASCENDING NODEIAM! - 
KEY4=1 FOR TIME CORRECTION LOOP, =0 FOR NO CORRECTION, 

KEY5=1 FOR BROKEN PLANE INPUT, =0 FOR STANDARD INPUT, 

ISEARCH=1 FOR PARAMETER RUN ON PARTICULAR LAUNCH-ARRIVAL DATE PAIR,0=NO 
PER1,ELP1,GEE1,XI l.DAYl = FIRST VALUES TO BE INPUT WHEN SEARCHING 
ON A PARTICULAR LAUNCH-ARR I VAL DATE PAIR, 

PER2, ELP2,GEE2,XI2, DAY2 = LAST VALUES. 

KPER,KELP,KG££ ,KXI , KOAY = INCREMENTAL VALUE S ( I NTEGER ONLY! 
TOEST*ESTIMATED TIME FROM DEORBIT TO LANDING, =0 YIELDS COMPUTED TIME, 

ALL LENGTHS IN KILOMETERS, ALL ANGLES IN DEGREES. 
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WRITE(6,650» 

IF(KEYl.Ne.O)GO TO 2 

WRITE(6,875I 

MRITE(6»876) 

2 CONTINUE 

IF(KEY5,NE.UG0 TO 21 
24 READ(5»BPOATA> 

21 CONTINUE 

CALL CALJULCEMJO.EF JD»WND»FD»EI 
CALL CALJUL(AWJD,AFJDfWNO.FD,XM» 

DO 700 I=ltILO,IJD 
WRIT£(6,500> 

00 700 J=1,IAD,JJ0 
DATeE=FLOAT(I-lI+EWJD+EFJO 
DATGM=FLOAT( J-D + AWJD+AFJO 

EWDI=IFIX(DATEt» 

EFCI=DATEe-EWDI +.0000001 
AWDI=IFIX(DATEM) 

AFDI=0ATEM-AW0I +.0000001 
CALL JULCAL(E1,EWD1,EF0I,0» 

CALL JULCAL (XMl .AWOI ,AFDI .0) 

1F(KEY5.£0.1)G0 TO 22 
IF(KEY1.NE.0»G0 TO 4 

CALL PLUG(DATEE.DAT=M,C3tDLA,VHP.0PA,RAP,KK) 

GO TO 5 

4 CALL PLUG1(0ATE£,DATEM,C3,DLA,YHP»DPA,RAP,TT.RAL,0C0M,PR0BPER,PR0B 
lAP,PR08INC,TAL.TAA,HELANGfGDA,6RA,ZAP,ETS,ZAE»ETEf ZACf ETC.KKI 

5 CONTINUE 

IF(KK.EO.O»GO TO 700 

IF(C3.GT.C3MAX)G0 TO 700 
IFCDLA.GT.OLAMAXIGO TO 700 
IF(VHP.GT.VHPMAX»GO TO 700 
GO TO 23 

22 CONTINUE 

23 CONTINUE 


IF (ISEARCH.NE.OIGO TO 12 
JPER=1. 

JELP=1. 

JGEE*1, 

JXI =1. 

J0AY=1. 

12 CONTINUE 

IFdSEARCH.NE.llGO TO 13 

JPER=ABS(PER1-PER2)+1 

JELP=ABS(ELP1-ELP2)+1 
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JGiE=ABS(GEEl-GEE2)+l 
JXI »ABS( XII- X12)+l 
JDAY=ABS(DAY1-DAY2) +1 

13 CONTINUE 
IF(KEY1.NE.2)G0 TO 15 

WRITE(6.900)E1(2) tEl(3) lEKl) t E 1 ( 4 ) .E 1 ( 5 ) » E 1 ( 6 > » XMl ( 2 ) t XMl ( 3 I , XMl ( 
lU ,XM1 (4) ,XM1(5) ,XM1(6I ,OATEE,DATEM 
15 CONTINUE 

DO 700 I1=1.JPER,KPER 
DO 700 12=1 .JELP. KELP 
00 700 I3=1,JGEE.KGEE 
DO 700 14=1, JXI ,KX1 
DO 700 15=1 ,JOAY,KDAY 

IF (ISEARCH.NE.UGO TO 14 
PER=PER1+I1-1. 

ELF=ELPl+I2-l. 

GEE=GEE1+I3-1. 

XI = XI1+I4-1. 

DAYS=DAY1+15-1. 

14 CONTINUE 

KEY2=1 

CALL COMPEL(PER,ELP ,GEE,XI ,CAPW,XITW,DATEM,HA,HP,AE,EE,DAYS,RS,KEY 
12,KEY3,EL0NP) 

EPRAD=HP + RS 
1F(KEY2.EQ.C)G0 TO 700 

CALL FUDGE! XJ2, XI ,AE,EE,OELCOM,DELSOM,PO,XITW,CAPW,OE,WE ,U,DAYS,RS 
1> 

IF(KEY1.NE.C >G0 TO 6 
M=10 

CALL 0EB00ST(AE,EE,X1 , WE , OE , Ft , U , DP A , RAP , VHP , RPM I N , I NO , AZ , EH , Z I , W 
1Z,0Z,FH,BT,BR,0ELTV ,TPERE ,M» 

GO TO 7 

6 M=1 

CALL DBST1(AE,EE,XI ,W£,OE,FE,U,DPA ,RAP ,VHP , RPM I N , 1 NC , AZ , EH , Z I , W 
1Z,0Z,FH,BT , BR,OELTV,TPERt,TPERH,a,HPRAD,DBRAD,VXH,VYH, VZHfVXE.VYE, 
2VZE,VXD,VYD,VZ0,0ECHP,RAHP,DECEP,RAEP,VDBH,VDBc,PLANE,M! 

7 CONTINUE 

♦♦♦♦♦♦transformation OF BT,BR FROM EQUATORIAL TO ECLIPTIC. 

DR=. 017453292519943 
XS=COS(OPA*OR )*COS (RAP+DRI 
YS=COS( OPA*DR ) *S IN ( RAP*DR > 

ZS=SIN(DPA*DR» 

CALL RECEO(DATEM,0. ,0 . , 1 . , XEQ, YEQ, ZEO) 

CALL REOMEO(DATEM,XEO,YEQ,Z£0,XMEQ, YMEQ,ZM£0,DM,RM) 

CALL CROSS (XS.YS.ZS ,XMEQ, YMEO, Z ME 0,TX , TY , T Z , PRODT I 
CALL CROSS(XS,YS,ZS,TX,TY,TZ,RX,RY,RZ, PRODR) 
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C 

c 

c 

c 


c 


c 


CALL CROSS ( XSfYS, IS fC. t 0. t 1 . » Tt X , TEY t f PRODET) 

CALL CROSS (XSfYStZS t TEX ,TE Y , T6Z ♦ RE X ,REY , RE Z f PRODER ) 

CALL D0T<TX*TYf TZfTEX,TEY,TEZ,ANGl) 

CALL DOT( TX»TYf TZf REX,REYrREZ,ANG2) 

CALL OOT(RX,RY,RZrTeX,TEY,TEZf ANG3I 
CALL OOT(RXtRY,RZ,REX ,REY » REZ ♦ A NG4 ) 

BDT=BT 

BDR=BR 

BT-0OT*COS( ANGl^OR) ♦BDR^COSf ANGZ^OR) 
BR=BOT*COS<ANG3*OR)+0OR«COS(ANG4tOR) 
end of TRANSFORMATION. 

VfcXC£SS=OELVMAX-OELTV 

IFCDELTY.GT.DELVMAXIGO TO 700 

CALL SEESE(OAT£MtPDf A£,EEf XI,WE,3Ef UfRS,EX,EY,EZf IS,Ic f TSUNINtDURS 
lUN,TASIN,TASOUT,TEARINtDUREAR,TAEIN,TAEOUTf DELSOM,DELCOM,OAYS) 

IF( (KEYl.NE.l).OR. ( KEY4.NE.U )GO TO 17 

CALL £VENTS(EL3NP,SPECL0N, AE»EEf FEf OAYS.PO,DATcM,PER,FPA,R£NTRY,RS 
l,rOEST,U*TAOEOR0f TIMLANDfTIMD£OR»TIMOBST,VDEORBf KK) 

IF(KK.EO.O)GO TO 700 
0BW=IFIX(TIM08ST) 

OOW=IFIX(TIMDEOR) 

XLW=IFIX(TIMLAND) 

08F=TIMD8ST-DBW 
OOF=TIMDEOR-OOW 
XLF=TIMLANO-XLW 
CALL JULCAL(DBST»CBW»DBF,0) 

CALL JULCAL(0£0R,DQWf OOFtO) 

CALL JULCAL(XLANOrXLWf XLFf 0) 

GO TO 20 
17 CONTINUE 

SPECLON=ELONP 

TACe0PB=0. 

VDEORB=0. 

20 CONTINUE 

IF(KEY1.E0.2)G0 TO 16 
IFCKEYl.EQ.l )GO TO 8 

WRITE (6,800 )E1( 1) »tl(2) ,E1 (3) ,XMU1 >,XM1(2 ) f XM1(3) , C3 , DL A , D£ LTV , FE 
1 ,FH,TPERE,IS, IE,TSUNIN,DURSUN,TASIN,TASOUT ,TEAR I N, DURE AR ,T AE IN , TAE 
20UT 
GO TO 9 
8 CONTINUE 

IFCKEY4.NE.UG0 TO 18 

WRITE! 6,899>E1( 2) ,E 1( 3) ,E 1 ( 1 ) , L 1 ( 4 ) , E 1 ( 5 ) , E 1 1 6 ) , XMl { 2 ) , XMl ( 3 ) , XMl C 
11) ,XM1 (4) ,XM1 (5) ,XM1( 6) ,DBSr(2) ,DBST(3) , DBS T ( 1 ) ,OB ST ( 4 ) , 08 ST ( 5 ) ,08 
2ST(6) ,DE0RC2) ,OEOR( 3) ,DE0K(1) , DEOR( 4) , DEOR( 5 ) ,0E0R(6) , XL AND C 2 I , XL A 
3ND(3) tXLANDd) , XL AN D( 4 ) , XL AND { 5 ) , XL AND ( 6 ) , DATE E , DATE M, T I MOBST , TI MD 
4E0R,TIMLAND 
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GO TO 19 

18 CONTINUE 

WRITEC6,900)E1(2»,E1(3» ,El( II ,El(4> ,E1 (5) ,£1(6) ,XMl(2 J f XMl (3> ,XM1( 
11) ,XM1(4),XM1(5),XM1(6) ,OATcE,DATEM 

19 CONTINUE 

WRITE(6,901 IOLA,RAL,C3,TT,OPA,RAP,VHP,GOA,6RA,OCOM,ZAP ,ETS ,IA=,ETE 
ItZAC.ETC, PROSPER, PROS AP,PR06INC,TAL,TAA,HeLANG,B,BT,BR 

WRITE(6,902)AZ,EH,ZI,OZ,WZ,FHI,TPERH,HPRAD,OECHP,RAHP, VDBH, VXH,VYH, 
IVZH 

WRlT£(6,903)A£,£E,KI,OE,WE,Fe,TP6RE,TADEORB,VD£ORB,EP^AD,OECEP,RAE 

1P,PD,V0BE,VXE,VY£,VZE 

WRITE(6,904IO£LTV,VEXCESS,CBRAD,PER,VX0,VyD,VZD, PLANE 

WRITE(6,905)IS,TSUNIN,0URSUN,rASIN,TAS0UT,IE,TEARIN,DJREAR,TAEIN,T 

lAEOUT 

WRITE(6,906)SP£CL0N ,ELP,GEt,OAYS 
GO TO 9 
CONTINUE 

WRITE (6,910 ) PER.ELP ,GEE,XI ,OAYS,OELTV 
CONTINUE 

CONTINUE 
CONTINUE 

IF(KEY5.E0.1 IGO TO 24 
GO TO 1 


F0RMAT(*0*) 

FORMAT (*l*l 

F0RMAT(1X,6F3.0,2F6.2,1X,F5.3,2F6. 1 , IX , F9 . 2 , 2 1 3, 2F 9. 2 , 2F8. 2 , 2F9 . 2 , 
12F8.2) 

875 FORMAT(* LAUNCH ARRIVAL C3 OLA DELTAV F DBST F DBST OBST T 

IIM SO EO TSUNIN OURSUN TASIN TASQUT TEARIN DUREAR TAEI 

IN TAEOUT*) 

876 FORMAT!* DATE DATE ELIPSE HYPERB FROM P 

lER*/) 

899 F0RMAT(20X,*LAUNCH OA IE* , 12X , *ARRI VAL DATE* ,12K , *DEBOOST TIME*,12X 
1, ♦DEORBIT TIME*, 12X, ♦LANDING TIME*,/,* CALENDAR* , 6X , 6F 3. 0 , 6X , 6F3. 0 
2,6X,6F3.0,6X,6F3.0,6X ,6F3.0,/,* JUL I AN* , 5X , FI 8 . 2 , 6X , FI 8 . 2 , 6X , F 1 8. 2 
3,6X,F13.2,6X,F18.2, /) 

900 F0RMAT(21X,*LAUNCH DATE*, 13X , *ARRI V AL DATt*,/,* C AL ENDAR*, 7X , 6F3, 0 
1,7X,6F3.0,/,* JULI AN*,7X,F 18. 2,7X,F18. 2, /) 

901 FORMAT!* INTERPLANETARY FLIGHT PARAMETERS*,// 5X,*DLA »*,cl6.8, 5X 

1, *RAL =*,E16.8, 5X,*C3 =*,E16.8, 5X,*TRIP TIME =*,E16.8,/ 5X,*AR£0 

2. DEC, S-VECTOR =*,E16.8, 5X,*ARE0. R.A. S-VECTOR =*,E16.8, 5X,*HY 
3PER, EXCESS VEL. =*,E16.8,/ 5X,*GE0. OfcC. S-VECTOR =*,£16.8, 7X,*G 
4E0. R.A, S-VECTOR =*,E16.8, 5X, *COMUNI CATI ON DIST. =*,E16,8,/ 5X,* 
5ZAP =*,£16.8, 5X,*ETS =*,E16.3, 5X,*ZAE =*,£16.8, 5X,*bTE =*,E16.8 
6,/5X,*ZAC =*,E16.8, 5X,*ETC =* , E 16. 8 , /5X , *PROB£ PERIHELION =*,E16. 
78, 10X,*PROBE APHELION =*,E16.3, 7X,*PR0BE INCLINATION =*,E16.8,/5X 
8,*LAUNCH TRUE ANOM. =*,E16.8, 5X,*ARRIVAL TRUE ANOM. =*,E16.8, 5X, 
9*HELI0, ANGLE TRAVEL =*, E16. 8 ,/ 5X, *B- VECTOR MAGNITUDE =*,E16,8,15X 


16 

9 

C 

600 

7P0 


500 

650 

800 
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$,♦6 DOT T =*»E16. 8» 17X»*B DOT R =*tE16.8»/» 

902 FORMAT** ELEMENTS AND DEBOOST PARAMETERS - HYPERBOLA*»// 5Xt*A =♦ ♦ 
1E16.8, 5X,*E =*,£16.8, 5X,*I =*,E16.8, 5X,*CAP. OMEGA =*,E16.8, 5X, 
2*0MEGA =*,E16,8,/ 5X,*DBST TRUE ANQM, =*,E16.8, 7X,*DTST TIME =*,E 

316.8, /5X,*PER. RADIUS =*,E16,8, 5X,*PER. DEC. =*,E16.8, 6X,*PER, R 
4. A. =*,E16.8,/5X,*V AT DBST =*,E16.8, 6X,*VX AT DBST =*,E16.8, 5X, 
5*VY AT DBST =*,E16.8, 5X,*V2 AT DBST =*,E16.8,/» 

903 FORMAT** ELEMENTS AND OEBOOST PARAMETERS - ELLIPSE*,// 5X,*A =*, 
1E16.8, 5X,*E =*,ei6.8, 5X,*I =*,E16.8, 5X , *CAP. OMEGA =*,E16.8, 5X, 
2*0MEGA =*,E16.8,/ 5X,*D8ST TRUE ANOM. =+,EI6.8, 7X,*DBST TIME =*,E 

316.8, 5X,*00RB.T. A. =*, £16.8 , 5X , *VD0RB =*, E 16. 8 ,/5X , 

1 *PER. RADIUS =*,E16.3, 5X,*PER. DEC. =*,E16.8, 6X,*PER. R 
4. A. =*,E16.8, 9X,*PERI0D =♦, El 6 . 8, /5X, *V AT DBST =*,E16.8, 6X,*VX 
$AT DBST =*,E16,8, 5X, 

5*VY AT DBST =*,E16.8, 5X,*VZ AT DBST =*,E16.8,/) 

904 FORMAT** DEBOQST MANEUVER PARAMETERS*,// 5X,*DELTA V =*,E16.8, 5X, 
1*EXCESS DELTA V =*,E16.8, 5X,*RADIUS =* , El 6 . 8 , 5X ,*PER =*,E16.8,/ 5 
2X,*VX =*,ei6.8,13X,*VY =* , E16. 8 , 17X, *VZ =*,E16.8, 9X,*PLAN£ CHANGE 
3 =*,E16.8,/) 

905 FORMAT** OCCULTATION PARAMETERS*,// 5X,*1ST ORBIT, SUV=*, 1 3, 8X,*TI 
IME, SUN =*,£16.8, 7 X, *OURATION , SUN =*,E16.8,/ 5X,*TRUE ANOM., SUN 

2 IN =*,E16.8, 7X,*TRUE ANOM., SUN OUT =*,E16.8,/ 5X,*1ST ORBIT, EA 

3RTH =*,I3 , 5X,*TIM£, EARTH =*,£16.8, 5X , *OURAT I ON , EARTH =*,E16 

4.8, / 5X,*TRUE ANOM., EARTH IN =*,E16.8, 5X,*TRUE ANOM., EARTH OUT 
5=*E16.8,/» 

906 FORMAT** LANDING POINT PARAME Tt RS* , //5 X ,*L3NGI TUDE =*,E16.8, 5X,*L 
LATITUDE =*,E16.8, 5X,*SUN LIGHTING ANGLE* , E 16. 8 , 5X,*DAYS =*,F3.0, 
2/*l*» 

910 FORMAT*/, IX, *PER =* ,£12.4, 3X,*ELP =*, E 1 2, 4, 3X , *GEE =*,E12.4,3X,*XI 
1 =*,E12.4,3X,*DAYS =* ,E1 2 ,4, 3X , *DELTA V =*,E12.4) 

END 
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SUBBOuriNfc PLUG(DATEEf DATEM,C3f DLA,VHP.DPA,RAP,KK> 

USUN=1. 3271411E+11 

CALL EfcARThCDArtt »X£ ,YE»ZE,OXE,DYE,OZE \ 

CALL fcMARS(CATEM,XM ,YM,ZM,DXM,DYM,OZM) 

TT=DATfcM-CATEE 

CALL LAM8RT (Xt » YE » Z £ , XN » YM, ZM, TT#24.*3600 . , A , E » XI , W » 0 , TAl i TA2 i USUN 
1 .KK) 

IF(KK.LQ,0)RETURN 

CALL CCNCAR (AtC,XItWtO,TAl,Xl,Yl. ZlfDXl,DYltDZlfUSUN) 

DXlt=DXl-DXE 

DY1£=CY1-DYE 

OZi£=CZl-DZE 

C3=DX1E=*=*2+DY1E**2 + CZ1E**2 

CALL RtC£0( CATcEfDXlE.DYltfCZlc ,X£Q,YtOiZEQ) 

CALL LArLNG(X£C,YEQ »ZEO,DLA,RAL) 

CALL CCNCAR (Af E.XI ,W,UfTA2fX2tY2t Z2tOX2tDY2tOZ2.USUN) 

DX2M=DX2-DXM 

DY2M=DY2-DYP 

DZ2R=CZ2-DZP 

VHP=SORT(DX2M**2+DY2M**2+DZ2M**2) 

CALL RLCEQI CATEM,CX2M,DY2M,OZ2M ,XEO»YtO. ZLO) 

CALL RcOMcOlDATEM, XEO.YEOf ZEQf XMEQ.YMEOt ZMECf DPA,RAP) 

RETURN 

END 


SUBRCUTINc PLUGl ( DA TEE t C ATcM, C3 i DL A , VHP f OP A .RAP , T T , RAL » OCOM , PROBPE 
IR.PROBAP.PROBINC, TAL, TAA.HELANG.GDA.GRA.ZAP.ETS , ZAE , ET E , ZAC , E TC , KK 
2 » 

USUN=1. 3271411E+11 

CALL EEARTh (CATEL, XE, YE.ZE.DXE.DYE.OZE ) 

CALL tMARS( CATEM ,XK ,Y M , ZM , OXM, DY M, 0 ZM » 

CALL LCARTH(DArEM, XEA.YEA .ZEA.DXEA.DYEA.OZEA) 

XMt=XtA-XM 

YMt=YEA-YM 

ZMc=ZfcA-ZM 

DCCM=SQRT ( XPE*XME+YME*YME+ZME*ZMc) 

TT=OATtM-OATtt 

CALL LAMBKr(XE,YE,ZE,XH,YM,ZM,TT<‘24.*3600.,A,£,XI,W,0,TAl,TA2,USUN 
l.KK) 

IFIKK.tQ.OJ RETURN 

PROBPlR=A-A*E 

PRQBAP=A+A*E 

PROaiNC=XI 

TAL=TA1 

TAA=TA2 

hcLANG=rAA-TAL 

CALL CCNCAR I A.E ,XI .W.G.TAl.Xl.Yl.Zl.DXl.DYl.DZl .USUNJ 
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CXlt=DXl-DXf: 

DYlfc = OYl-DY£; 

DZlc=CZl-DZb 

C3 = CXlfc=f‘*2 + DYlE**2+DZlE**2 

CALL RfcCCO( CAT££,DXlE,DYlt.DZll; ,XEO , YEO, ZEQ ) 

CALL LArLNG(XfcC.YfcCfZ£QrDLA,RAL) 

CALL CCiMCAR{A,E»X I , W t Q ,TA2 » X2, Y 2 1 Z2 1 0X2 »0Y2 ,DZ 2 . USUN » 

0X2M=DX2-DX|v 

DY2M=OY2-OVM 

0Z2N=DZ2-CZM 

CALL OCT {-XNi,-YM ,-Zl^iOX2M,DY2M,OZ2M,ZAPl ) 

CALL DOT( XMe,Yy£,ZMt.DX2M,DY2M,0Z2M,ZAEl> 

VHP = SORT ( DX2M<=*2 + DY 2M«*2+DZ2M**2 ) 

CALL RuC^O( CAT cM,DX 2M,DY2M,DZ2M ,XtO,YfcO»ZCQ) 

CALL LArLNG( XEC.YcC ,ZEg,GDA,GRA> 

CALL RtOMfcO (DATEM,X£OtYfcO,ZEO,XMcOt YML-Q»ZMEO,DPAfRAP) 

SX=XMLC/VhP 
SY=YMEC/VHP 
SZ=ZMLQ/VhP 
SRT = SORI (SX'i'SX + SY’i'SY) 

TX=SY /SRT 
TY=-SX/SRT 
TZ = C. 

SX1=0X2M/VHF 
SYl=DY2M/VhP 
SZ 1~0Z2M/VHP 

SRT1=SQRT (SX1*£X1+SY1*SY1 ) 

TX1=SY1 /SRTl 
TYl=-SXi/SPTl 
TZ1=0. 

CALL RcCEOi CATfcM,TXl,TYl,TZl.TXEQ,TYEQiTZfcQ) 

CALL RLCML0(CAT£M,TXEg,TYEQ,TZEQ,TXME0iTYME0iT2MEQf0ECf RA) 

CALL DCT ( rx,TY,TZ,T XMEQtT YMEQ,TZM£Q, ERROR) 

CALL CROSS ( SXf SY.SZ .T X , T Y . TZ » RX » RY , RZ t RMAG ) 

CALL VtCTOR (DATLMi XI, X2f X3,X,X, XiSUNX,SUNYiSUNZ»fcXiEY,EZ,CX,CY,CZf 
lA) 

SLNS=SX*SUNX+SY*SUNY+SZ*SUNZ 

SUNr=TX#SUNX+TY=i'SUNY+ rZ*SUNZ 

SUNR=RX*SUNX+RY*SUNY+RZ*SJNZ 

EAS sSX’i'LX +SY>l<tY +SZ*EZ 

LAT =TX*aX +TY*2Y +TZ*tZ 

LAR sRX’S'tX +RY*cY +RZ*CZ 

CAS =SX*CX +SY*CY +SZ*Cl 

CAT =TX#CX +TY*CY +TZ*CZ 

CAR =RX*CX +RY*CY +RZ*CZ 

CALL LATLNG (SUNT ,SUNR, SUNS, SOEC, SRA) 

CALL LATLNG(£AT ,EAP ,E AS ,E DEC , £R A ) 

CALL LATLNG (CAT ,C AR ,C AS.COEC ,CRA ) 
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ETS=SRA+180. 

ETE=EKA+18C, 

tTC=CRA+180. 

ZAP =90^-SDfcC 

ZAE=9C.-£DLC 

ZAC=9G.-CD£C 

IF(ABS(ZAPl-ZAP).GT.l.)WRiTE(6,100)ZAPl,ZAP 
IF (A8S< ZAEl-ZAE).GT.l. ) WRITE (6,200 > ZAhltZAE 
100 FORMAT (2E16. 8,* ERROR IN ZAP«) 

200 F0R^AT(2Z16,8,* ERROR IN ZAE^) 

RETURN 

END 


SUBROUTINE COMPEL C PER , ELP , GEE , X I , CA PW , X I TW , DAT EM ,HA , HP , AE , EE , 0 A YS 
1RS,KEY2,K£Y3,ELCNP) 

ALL ANGLES INPUT IN DEGREES AND OUTPUT IN DEGREES^ 
ANGLb(X)=AMOD(X,3 60.)+180.-SIGN( 180. ,X) 

XSIN(X)=S1N(DR*X) 

XCCS(X)=COS(DR«X) 

IF(ABS(XI ).LT.ABS(ELP) )G0 TO 50 
DR=. 017453292519943 
C JUL=CATcM+EAYS 
1C0D£=4 

CALL VECTOR ( D J UL , E L S , E LG NS , DEC E , R AE , DLCC , R AC , S X , S Y , S Z , E X ,£ Y , E Z , CX 
lCY,CZ,ICODc) 

IF (ABS(SlGN(GEt ,LLP )+ELS) .LE.ABS(ELP) )G0 TO 50 
RD=57.295779513C82321 

ARG=( XCQS(GbE )-XSIN(ELS»=«'XSIN(ELPl ) / ( XCOS ( EL S ) *XCOS ( E L P ) ) 

IF(ABS( ARG) .GT.l. )GC TO 50 
ARC1=ASIN( XSIN(LLP)/XSIN(Xn )*RD 

ARC2 = ASIN( ( XSINC tLP)/XCOS (ELP) ) « ( XCOS ( X I) / X S I N ( X I ) ) )*^D 
GO TO ( 10, 20,30, 40)KEY3 
10 CONTINUE 

ELQNP=( ACOS (ARG) ) *RC+ELCNS 
XITW=PER+180.-ARC1 
CAP^=EL0NP-18C. +ARC2 
GO TO 100 
20 CONTINUE 

tLONP=( ACOS(ARG ) )=«'RD+ELGNS 
XI TW=PlR+ARC 1 
CAFw=ELGNP-ARC2 
GO TO 100 
30 CONTINUE 

ELONP=ELONS-( ACCS (ARG) )’«'R0 
XITW=PEK+18C.-ARC1 
CAPW=EL0NP-180.+ARC2 
GO TO 100 
40 CONTINUE 
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fcLONP = ELONS-MCCS(ARG) )*RU 
XI TW=P^R+ARC1 
CAPW=t:L0NP-ARC2 
60 TO 100 
ICO CONTINUE 

XlTh=ANGLE(XITWl 
CAPW=ANGLfc (CAPW) 

PA=RS +HA 
RP=RS +HP 
At= (RA+RPI/2. 

££ = (RA-RP»/(2.*Ai:) 

RETURN 
50 Kt¥2=0 
RETURN 

e:no 


SURROUriNE FUCGK XJ2f XI ,A£ ,tb,DELCOM,DELSOM,PO,XITW,CAPWtOE, WEfUtO 
lAYS.RS) 

XN=S0RT(U/A£**3) 

PI = J. 1^1 5<i2£ 53 58975 3 
DK=PI/18u. 

C1 = (RS / ( Ac*( l.-£E*cE m*«‘2 

DlLSC«=6.*PI*XJ2*C1«( 1.-1.25*SIN(XI*DR )**2) 
0CLC0M=-3.*PI*XJ2*C1*C0S( XI*DR) 

cNBAR=XN*(l. + 1.5*Cl*XJ2*SQRT( 1 . -EE *EE >« ( 1 . - 1. 5*S IN ( X I *DR ) **2 ) ) 
PD=2.*PI /ENSAR/8640Q. 

B=CAYS/PD 

Ol;=CAPW-( B*CLLCCMJ*180,/PI 
Rb=Xl TW-(B*0ELS0M)*180./PI 
RETURN 
END 


SUHROUTINL .SElSL(CATEM»PDfAi:»£E»XI»W£»OE»U»RS»cXfEY»EZ»IS»IE»TSUNI 
lN,OUKSUiM, TASIN, TASOUT,TcARlN,DUREAR, TAEINtTAEOUr ,0£LSQM,DELC 
2UM,DAYS» 

IS=C 


lt=0 

TSUNIN=j. 

DURSUN=0* 

TASIN=3, 

TAS0UT=0. 


TE ARIN=0. 

DUKcAR=0. 

TAcIN=0* 

TAt0UT=0. 

ISTGP=CAYS/PD+1 
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DO 20 J=ltISTCP 
TlM£=CATtM+PD*FLOAT ( J-1) 

CALL VfcCTORdIMC ,DECS tR AS t DfcCh » RAE , DLCC t RAC » SX , S Y t SZ , tf X , £Y , L Z t CX , C 
lYtCZ,4) 

CALL OCCULT (A£ » E£ t X 1 » WL + PD*FL JAT ( J-1 ) ^DLLSOM ,0t: + PD«FL3 AT ( J-1) «DELC 
10MtU,RS» SX,SY,SZf UUKSUfTSUNl f ALTl tTASi » DEC 1 • R A 1 , T2 , AL T 2 » TASOU , D£C2 
l»RA2f KS) 

CALL OCCULT (A£ .££ t X 1 , WE + PD*FLOAT ( J-1 ) ^DELSDM, Oc+PO*FLDA T< J- 1 KD £LC 
ICMfUf RS*EX,ey,LZtOURL,TLARI f ALTlf TAdI , DLC 1 , P A 1 , T2 1 AL T2 r FAcCUf DEC2, 
lRA2»Kc ) 

IF( IS.NE.O)GO TC 10 
lF(KS.fc0.1)G0 TO 8 
GO TO 10 
8 IS = J 

TSUMN=TSUNI 

DURSUN=OURSU 

IASIN=TASI 

TASCUT=TASOU 

10 IF( It . Nc. 0 )GO TO 15 
lF(Kc.E0.1)G0 rC 11 
GO TO 20 

11 IE = J 
rLARIN=TLAkl 
DURfcAP=DURE 
rAcIN=rAcl 


TAECUT=TAtOU 

15 IFdc.NE.O .AND. I S .NE . 0 ) RETURN 
20 CCINiTINUt 


RETURN 

END 


SUE ROUT I Nt cVcNTS( irLCNP.SPtCLUNt Act EL » FE r OA YS , PD , OAT EM t PER t FPA t RE N 
1 TRYfRSt TDlST ,U»TADL 0R3,TIMLAND, TlMDtORr TIMOBST f VDhOKB, KK ) 

COMPUr.-S EVENT TIMES FCR LAND I NG t DEOK B I T ANG DcBOOST, GIVi^N LANDING 
PUINT AND cLLM£^TS OF ELLIPSE AND TRUt, ANOMALY OF DtORBIT. 

currl-Ct fur time of day. 

ANGLE ( O = AMU0( X, 360. I +1 80 .- S I GN ( 1 80 . r X ) 

PMDCT=350. 891962 

HA= 145. 845 + 3 50. 89 1 9 62* I DA T EM + DA YS- 24 1 8 322 . ) 

DLLrLOK=LLGNP-SP£ CLCN-HA 
DELTLUN=ANGLe( DELTLGN) 

CcLTJG=DcLTLON/PMOCT 
TIMLAND = DA T tM+C AYS+ CfcLT JD 


IF(TDESr.NL.O.)GG TO 1 
Pi.RG=-PEK 

CALL CCNFPAI AE ,EL*TADEORB,PERO»RS»RENTRY,FPA,Ut AL,EL,FLOfFLO»VDLOR 
IBtThfcfA.KK) 
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IF (KK.iiO.OI PETjRi'J 

CALL TCGNIC (U»l;L» AL »FLC.TDEOR) 

CALL TCGNlCIUtcLf AL .FLC.TLANOJ 
TDH0k=T9t0k/86ACC. 

TLAM)=TLAND/864CG. 

PDT = 2.*3,lA15C2653b=»S0kn AL*AL*AL )/SOkT(U) 
PL)T=PDT/864C0. 

IF ( IDLUK.GT.C. ) CLLDt:CR=7LANO+PDT-TDEOk 
IF ( TOi;Uk.LT.:j. ) CLLDcCR=TLAND-TDt:OR 
GU TO 2 

1 CCNTINbc 
DLLDi.Ok=TDlST 

2 CCNTINUfc 

f 1 MDcOk=TI MLANG-D tLOLOk 

CALCULATE DLBCGST TIMc, 

6ACKUP=DAYS/PD 

Df.LTIk=IFIX(BACKUP) 

CALL ICGNIC(U,2i.,AL »Fu,rPtRS) 

CALL TCUiNjlC (U,tf ,A£tTADEJKB,TDL-OR3> 

TPhRfc= rPtRh/ecACC. 

TU:=CKB= TDL0kB/8 64UC . 

IF ( rut.GR8.GT.C. ) D2 L CD B= IDt ORB- T PERt 
IF ( ID^URB.LT.v'. ) D2LCCB=TDcOR3-TPtkL- + PD 
TlMUBSr =ri FD2Gk-U2LTIk-DL.LQDB 
R£ TUkINi 
lND 


SUB ROUT INI LLBOGSn A1 »L1, 1 1 , W1 , 01 , F 1 , U , L AT S . LONS , VI NF , R PMI N , INC.AZ 
$,t Z ,IZ,NZ,UZ,FZ, BT , BR,DLLTV»TPLRd,M » 

REAL Il,IZ,LArS,LCNS,NX,NY,NZ,N,IZP,lZM 

CATA DR ,kU, PI/. 17 4532 S25 1994 3F-1, 57.295779513082321, 3. 141592653589 
4793/ 

ANGLc (X)=AMCC(X,2.=>P1 ) +P I - S I GM P I , X ) 

Cl RCNSICN L)V(36C ) ,T A( 3bO) ,HYP( 363,6) 

OIPCNSILN PX(3) ,PY( 3,6) 

CLAT=CuS(OK=!‘LATS ) 

SLAT=SIN(Oh=i'LATS) 

CLCN=CCS(DR*LCNS) 

SL0N=SIN(DR«L0NS) 

SX = CLAT*CU,N 
SY=CLAT*SLCN 
SZ=SLAT 

DO 1 I=l,3bv.,M 
TA(1)=FL0AT (1 )-18C. 

F = Ck<‘TA( I ) 

CWF=CCS ( DR»kl+F ) 

SNF=S1N(0K*R1+F) 

CI=COS(DR*I 1 ) 

S1=S1N(0RX'11 ) 
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CU=C0Sl0R*01 ) 

SO = SIN(0R=!'Ol » 

RX=CwR*CO-SWF*Sa*CI 

RY=CWF*SU+SWF»CO*CI 

RZ=SWF*SI 

AZ=-U/VINF=i‘*2 

R0 = A1>!=( ) / ( l. + £l#COS (F) ) 

RS = RX>i‘SX+RY*SY + R2*SZ 
A=AZ*>i‘2 

B=R0**2=i'RS**2 + 2.*R0*A Z*RS-R0**2-2 . *AZ**2 +2.*AZ*R0 

C=2,=*=R0**2#RS-2.*RC*AZ*RS + 2,*R0**2 + AZ#*2-2, *AZ*R>j 

0V{I»=1,&20 

TfcST=B=f'B-4.*A*C 

1F( TfcST.LT.O. )GC TO 1 

OISC=SOKr(TtST» 

£ZP=SQRn (-e + DISC)/ (2.*A) ) 
fcZR = SORn (-B-OISC) / (2.*A) > 

IF(LZM,Lt.l.)i:2R=£:2P 
PhiP=ACOS( 1 ./£ZP) 

PH1M=AC0S( U/fcZR) 

FZP=ACQS( (AZ#( l.-LZP**2 )- RC » / { £ ZP*RO » ) 

FZR = ACOS( (AZ*( l.“tZR**2)“R0)/(L;ZM*R0) » 

IF(A8S(C0S(ANGLL( PHIP-FZP ) )-RS ) . GT . 1. 0-7 I F Z P=-F Z P 

IF (ABS(CQS( ANGLt (PMIM-FZM) )-RS) . G T . 1 . c-7 > F ZM=- FZM 

NX=RY*SZ-RZ*SY 

NY=RZ>!'SX-KX*SZ 

NZ = R X’!'SY-RY*SX 

N = SORT (NX*=i'2 + NY#*2 + NZ*#2» 

IZP=ACOS(NZ/N) 

IF ( ANGLc( P hIP-FZP J . GT .PI) IZP=ACQS (-NZ/N) 

IZR=ACOS( NZ/N) 

IF(ANGLt(PhIM-FZM) .GT.PI) I ZM=ACQS (-NZ/N ) 

IF( (IZP.Lti,Pl/2..ANC.INC .EQ.l) .UR. ( 1 Z P. G f . P 1/ 2. . AND. INC .E0.2) )2, 
S3 

2 LZ = i;ZP 
IZ=R0*1ZP 
FZ=RD*FZP 
FhI=RD>l'PHlP 
GO TO A 

3 EZ=EZM 
IZ=RD«IZN 
FZ=RD*FZM 
PHI=RO=('PHIM 

4 RPZ = AZ<'( l.-LZ) 

OV( 1 ) = 1.£20 

IF (RPZ.LT.RFMN )G0 TO 1 
WS=ASIN(SZ/SIN(DR*I Z) ) 

WZ=RD*WS-PhI 

IF (ABS(RZ-SIN(DR*(WZ+FZ) )>i'SIN(DR*IZ) ) . GT . 1 . t-7 ) WZ = 183 . -R D«WS- PH I 
DfcT=C0S (DR* IWZ + PHI ) ) * *2+COS ( DR* I Z ) **2*S I N ( DR* ( WZ + PH I ) )**2 
CC=(CCS(DR* (WZ+PHI ) )*SX+COS(OR*IZ)*SIN(DR*(WZ+PHl ) )*SY)/OET 
SO= (-COS (DR*I Z )*SIN (DR*(WZ+PHI ) ) *SX+COS(DR* (WZ*PHI ) ) *SY )/DET 
OZ=RD*ATAN2 (SO.CO) 
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hYP( I , 1 )=AZ 
rlYPd ,2)=i.Z 
hYPd »3) =1Z 
hYPd ,4)=WZ 
hYPd .5.)=OZ 
F1=PD*F 
hYPd ,6)=FZ 

CALL CCNCAPJAZ, LZ.IZfWZ,OZ.FZ.X,Y,ZfOX,OY,DZfU» 

CALL LCNCAP { Al, £1, Il,Wl,ai,Fl,X,Y,Z,VX,VY,VZ,U> 

CV ( n = SQh r { ( 3X-VX)*«2+( DY-VY)*«2+(DZ-VZ»#=S'2 ) 

1 CCMINU& 

1MIN=.'3 

ULLrv=l .£20 

DO 3 I = 1 , 36 3,M 

IF (DV(I I.Gf.DLLTViGC TO 8 

1MIN= I 

Ll-.LTV = DV( I ) 

8 CLMINLc 
1M1NM=1MIN-P 
IMIi\P=Ii>'.I.Y+K 

lF(IMII\h.LL-.0.0P.I?^INP.GL.361»G0 TO 6 

IF (DV d '^IN-h ) ,tC.l .t20.0K.DVdMIN + M ) ,£Q,1.E20 »GO TO 6 
PX ( 1 ) =TA( IMIN-M ) 

PX (2> =TA( IMN ) 

P X ( 3 ) = T A ( I K I N + M ) 

PYd .1 )=0V( IPIK-MI 
PY(2.1 )=0V( IMI!M> 

PY (3.1 )=DV( IPIh + P) 

CO 3 1=2,6 

PYd, I )=rlYP{IMN-P, I ) 

FY (2, I ) =riYP( ININ , I ) 

5 PY( 3, I )=nYP( IMN + N , I ) 

CALL PA KINIFl, O l LTV ,PX,PY( 1,11,0) 

LALL PAPIM(Fl,tZ,PX,PY( 1,21 ,1) 

CALL FAFIN(F1 ,lZ,PX,PYd,3) ,1) 

CALL PAPIN(F1,WZ,PX,PY(1,4),1) 

CALL PARIMF1,CZ,PX,PY(1,5) ,1) 

CALL PARIN( FI ,FZ ,PX ,PY( 1, 6) , II 
GO ri) 9 

6 Fl = lA ( m.'J ) 

AZ = hYP( IMIN.l I 
c Z = hYP ( I I'll K ,2 ) 

1Z=HYP(IMIN,3) 

WZ=hYP( IKIN,4) 

CZ=HYP( IMIN,5) 

FZ=hYPdMIN,6) 

9 B=-AZ*SCRT (eZ*f Z-l. ) 

Br=8^CCS(CP#IZ)/CLAr 
BR=8*SlNi(0H»IZ)*CCS ( DR* t LONS-QZ ) ) 

CALL TCLNIC (U,L1, A1 ,Fl,TP£Rr.) 

RLTURN 

tND 
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SUEROUriNc LBSri(Al,£;l,H,Wl,01tFl»U,LATS,L0NSFVINF,RPMIN,INC.AZ,E 
lZ.IZ»wZ,OZtFZ,BT,BR »OELTV t TPERb » TPBRH f B , HPR AO » DBR AO » VXri t VYH t V ZH t VX 
21-. . VYfe .VZs!: ,VXD»VYDf V ZDfDECHPf RAHP,DfcCrP .RAuP . VOBHrVDBt f PLANc.M) 

RfcAL 1 lfIZfLATS,LONStNX,NYtNZ,NiIZP,lZM 

DATA OR ,RD, PI /.I 7 45 32 ;25 19 94 3H-1, 57,29 5779513082 321 f 3. 14159265 3 589 
<793/ 

ANGLE (X)=AMCD (X ,2.*PI ) +P I- S IGN ( P I f X > 

OINEnSI JN DV(360) fTA(360) fHYP(36Df6) 

OIl^cNSION PX(3) fPY( 3f6) 

CLAI = CCS(CR*LATS ) 

SLAr=SIN(OR*LATS) 

CLCN=COS (DR*LCNS ) 

SLCN=SIN(DR#LCNS) 

SX=CLAT*CL!JN 

SY=CLAT*SLON 

SZ=SLAT 

CC 1 I=1f36GfM 
TA( I)=FLOAT(l)-ieO. 

F=CR«IA( I ) 

CVJF=COS(DR*m+F) 

SwF=SlN(OR*kl+F) 

C1=CQS(DR*11 > 

SI=SIN (OR*I 1 ) 

CO = COS(OR’i‘Gl ) 

SJ=S1N(0R»=0H 
RX = CWF*CJ-SRF*SO=«CI 
RY = CWF#SO + SWF’i'CO*C I 
R Z=SWF<‘SI 
AZ=-U/V1NF>!=*2 

R0 = A1^( 1.“L1=!'L1)/(1. + E1*C0S(FJ 1 
RS=RX«SX+RY*SY+RZ*SZ 
A = A Z**2 

B = R0**2l'RS=«'*2 + 2. *KG*A Z*RS-R0**2-2 . *AZ>!'’i‘2 +2,*AZ4 kO 

C=2, «R0**2*PS“ 2. *RC=»AZ»!<RS + 2.#R0**2+AZ**2-2. *AZ=!'RC 

ovm=i.t ;20 
TtSr = B*8-4.=»A>l'C 
lF<TfcST.Lr.C. )GC rc 1 
CISC=SOKr(TtST ) 
lZF=SCHr((-e+ClSC)/(2.*A)) 
tZM=SQRr( (-e-DISC) / (2,=«'A) ) 

IF(lZM.L.j.1.)LZR = 0ZP 
PmIP=ACOS( l./rZP) 

PHIN=AC0S( l./LZMI 

FZP = AC0S( (AZ*(1.-LZP**2)-RC )/(EZP*R0J ) 

FZR = ACOS( (AZ*(l,-tZW*=<'2)-RO)/(fcZM*RO» ) 

IF (ABS(COS( ANGLfc (PhIP-FZP ) )-RS) , GT . 1 . fc-7 > F ZP=-FZ P 

IF (ABS( COS ( ANGLE ( P4IM-FZM ) )-RS ) , GT , 1 , 8-7 I F ZN=-F ZM 

NX=RY4SZ-RZ*SY 

NY=RZ*SX-RX*SZ 

NZ = RX«SY-RY>!'SX 

N-=SCRT{NX*'!'2 + NY**2 + NZ*"t<2J 

I ZP=ACQS (NZ/N) 

IF( ANGLO ( PUP-FZP ) .GT.PI) IZP=ACOS (-NZ/N) 

IZR=ACCS( NZ/N) 

IF( ANGLE! Phi M-FZi^) . GT ,P I ) I ZM=ACOS {-NZ/N ) 
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IF ( ( I ZP. Ll. PI/2..AKC, INC .£Q.l ) . OR, ( I Z P. GT . P I / 2, . AND. I NC .£0.2112. 
43 

2 LZ=clP 

1Z = KU=!'1ZP 
FZ=kU*FZP 
PHI=KD*PniP 
GO TO 4 

3 hZ=uZM 

I Z=P0*1 ZM 
FZ = iW>!‘F ZM 
PHl = RU>!‘PnIM 

4 PPZ=4Z4<1.-CZ) 

DV ( I ) = 1 .c20 

IF (KPZ. L r.K PMIN )GlJ TO 1 
^S = /iSlNiSZ/SlMUk*l n ) 

WZ=PU*^S-Plil 

IF ( A>3i(RZ-S IN (L}P.’>= ( 1-1 Z + FZ) )*SIN(UR*IZ I ) . G T. I . fc-7) WZ = 180 . -RO*WS-PH I 

Cl T=CCS (Dk* (WZ + PmI ) )>i‘»2+COS (OR*IZ) #«2*SlN(DR#(.-iZ + PHl ) )**2 

cc= (ccs(Dk* ( wz+Prti ) t^i'Sx+cosoR^s^izi^siN (DR*( wz+PHi n*sy )/Dcr 

SC= (“COS (CP^I Z )*S IN (DR*(WZ+PHl I ) *SX+C0S ( DR# ( WZ+PH I n*SYI /D£T 

DZ=KD#ArAN2(SC.CC) 

hY P ( I , 1 ) = A Z 

FYP ( I ,2 I =: Z 

e Y P { I . 3 ) = I Z 

t- YP ( I ,4) =WZ 

I- YP( I ,5)=0Z 

F1=RD#F 

( YPd .o)=FZ 

CALL CCNCAF (AZ.tZ.lZ.WZ.DZ.FZ.X.Y.Z.OX.OY.OZ.U! 

CALL CLi\iCAP(Al,Ll.Ii.Wl,Ul,Fl.X,Y,Z.VX,VY,VZ,U» 

0Vm=S0KT( (DX-VX)**2 + ( DY- VY ) *#2+ ( D Z- YZ I ##2 ) 

I CCMINO.- 
IMIN=: 

Ct LTV = 1 .L20 

DO b 1 = 1 . 36C . M 

1F(0V( I ).f7r.D.LTV)GC TO 8 

1 M i \= I 

D^L rv = cv( I ) 
a CCM1NU-. 

I M i N K = 1 F' I N“ P 
I P 1 N P = I M 1 N + P 

IF (IMINM.LL ,C.0R.IMNP.Gi:-.361IG0 TO 6 
IKi;V(lMIN-P).l-y,l,i: 20 .GK.DV(IMIN+M).£; 0 . 1 .E 2 C)G 0 TO 6 
PX ( 1 ) =TA( IN IN-N ) 

FX (2 ) =TA( iMIN ) 

PX( 3)=TA( IMIN + M ) 

PY( l,l)=DVl iPI.N-M) 

PY{2.1)=DV( IMIN) 

FY(3,1 )=CV( IMN + MI 
DO 5 1=2.0 

PY( 1, I )=hYP( iPIN-M, I ) 

PY(2,I)=hYP(lMN,I ) 
iJ PY ( 3. I ) = hYP ( IMIN + M , I ) 

CALL PA kINIFI.OlLIV.PX.PY (1,11,0) 

CALL PAP IN( FI ,i.Z, PX ,P Y( 1 , 2 ) , 1) 
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CALL PARIfVHFl,lZ,PX,PY(l»3» ,11 
CALL PARIN(Fl,WZ,PX,py(l,4> ,1) 

CALL PAR1N(F1,GZ,PX,PY(1,5»,1> 

CALL PARIN(F1,FZ,PX,PY(1,6),1) 

GO ro 9 

6 F1 = TA( IMIN » 

AZ=HYPUMIN,1I 
EZ=hYP(IMIN,2l 
IZ=hYPi IMIN,3» 

WZ=nYP( 

OZ=hYP(IMlN,5) 

FZ=FYP( IMIN, 6) 

9 B=-AZ*SORr(cZ#eZ-l. ) 

BT = 8#CCS(0R«=IZ>/CLAT 
BR=8*SIN(DR*I Z»#COS (CR*(LQNS-OZ) ) 

CALL rCONIC(U,El,Al ,F1,TPERE) 

CALL TCCNIC(U,fcZ,AZ,FZ,TP£:RH) 

HPRAD =AZ-AZ*fcZ 

CBRAD=(A1-A1*L1#L1» /( 1. +fc l*COS ( DR*F 1 ) > 

CALL CCNCAH(AZ,EZ, I Z,WZ,OZ,FZ,X,Y,Z,VXH,VYH,VZH,U> 
CALL CC(NCAR(Ai,tl, Il,Wl,Gl,Fl,X,Y,Z,VXi:,VYt ,VZC,U» 
VDBH=SORT( VXH*VXh+VYH*VYH+VZH*VZri> 

VCBc=SORT (VXL*VXt+VYfc=i‘VYt. + VZ£*VZi:) 

VXL)= VXE-VXH 
VYC=VYL-VYH 
VZC=VZL-VZn 

CALL CCNCARI AZ, tZ, I Z,WZ,0Z,0., XPH, YPH,ZPH,DX,DY,[)Z,U) 
CALL CCNCAR(Al,51,Il,Wl,01,0.,XP£:,YPt,ZPt,DX,0Y,DZ,U> 
CALL LAFLNGI XPH ,YPH ,ZPH ,DlCHP,RAHP) 

CALL LArLNG(XP£,YPg,ZPc,DLCdP,RAi:P) 

WX£=SIN(DR*I1 I*S1N( CR*01> 

WYL=-C0S(CP»01 )*SIN(CR*11 » 

RZ£ = COS(OR*m 
UiXh=SIN(OR*IZ)*SIN( CR*OZI 
VvYh=-COS(OR*OZ)*SIN(OR*IZ) 

WZh=COS(DR*IZ> 

CALL OCr(UXt,fcY£ , WZ E , WXF , WYh , kJ ZH , PLANE ) 

RETURN 

lND 


SLERUUriN_ CAL JULCW JD,FJU,WNO,FD,X) 
DIMcNSlUN X(6» ,A( 12 ) 

050=2433282. 

YD = xm-48. 

YL=Y0/4. 

KYL=YL 

CK=KYL 

IF(YL-CK)1,1,3 
1 IF( X( 21-2. »4,4,3 

3 DS=CK 

GO ro 5 

4 CS=CK-1. 

5 DS=DS+365. *< YD-2. ) 
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DO 6 I = l » U 

6 A { I ) = i . J 
K=X(2 ) 

DO 7 1=K,12 

7 A { 11=0.0 

CS = DS + 31.*< A ( 1) +A ( 3 ) + a (5» + A (7 » + A(8»+A( 10 )+A( 12) » 
l + iC,>i=(A(<+)+A(6)+A(9 )+A(ll n+28.*A(2) 

CS=CS+X(3 1-1. 

WMC=DS 

FD=X( A) /2^.+X( 5 )/lA40,+X( 61/86400. 

IF (FD-.4999';99)9.6»a 

8 FJC=FD-.5 
WJD = 1 . 

GO ro lu 

9 FJD=FC+.5 

JD= ‘i , 

10 WJC=05^ +WJU+WN0 
HL rUKN 
LNC 


SUfiftOUTIN.'; CCjNCARI A fC.XI,W,0,F,X,Y,ZtDX,DY»DZ,lJ> 
DATA DR/. ji74632925i9943/ 

F R = DR F 
WFR = UR*i'( rt4-F } 

OR=CK*C 
XI fi=UK*X I 
Di-N=l.+^=i»CGS(FR ) 

{ 1. — )/Du.N 
V=SDRnU=i'(2./R-l./A) 1 
GAR = A TAN! ( FR )/DLN ) 

lrtFj8 = WFr<'»GAR 

CWF = COS ( WFR 1 
SwF=SIMWFk I 
SJ=SIy(Ok ) 

Cj = CljS (OR) 

S I = S I i'J ( X I R ) 

Ll=COS{XIk ) 

SwFG=SIN( aFGk) 

CWFG=CGS { wFGR ) 

X = R+(CaF«CC-SWF*SC<‘CI) 

Y = K’1‘(CwF*S0 + SWF’(‘CC*CI ) 

Z = K=i'SwF<‘S I 

CX = V>M-i»^FG*CC-CWFG*SQ*CI ) 

QY = V* (“Si'lFG*SO + CWFG*Ca’l<CI ) 

CZ=V*C'flFG*SI 

FE TURN 

LND 
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SOHRGUriNi CQNFPAt AC,tOtFO,PtRD,RSfREf FPAd,U»AL,eL,FLDtFLa.OtLVfTH 
lETA.KKI 
DK'.ENSIJN P(2» 

ANGLc (X) = AMCD( X, 3bC . » + ieO.-SIGN( 180, f X) 

DELV=0. 

DR=. 017453292519 
RD=d7. 2957795130 
KK=0 

ANG12 = ANGLL <PL-KC-F0 I 
S12=S1N( ANG12*0P) 

C12 = COS( ANG12*DR ) 

CFPA=CCS( FPA6*DF ) 

SFC=SIN( F0*CR) 

CFC=COS{FO>(=CR) 

Rl=AO«( l.-tC*cO )/( 1 ,+tO*CFO) 

R2 = RS 

V0=S0RT(U*(2,/R1-1./A0) ) 

FPAC=ASIN(EC*SFC/SCHT ( 1 , +2 .*1:0*0 FG +fcO*t:a ) ) *RU 
C 

A=-R2*R2-R1*R1+2. *R l*R2*C12+(Kl*R2*S12/Rfc/CFPA)**2 

B = 2.*(R1*R2*R2 + F1*P 1*R2-R1*R1*R2*C12-R 1*R2*R2*C12- (R 1*R2*S12 » **2/R 
lU 

C=Rl*Pl*R2*P2*(-2.+2.*C12+S12*S12) 

C 

CALL QAGRAT ( A ,S,C,P a » .P( 2 I ,KK » 

IFIKK.L-O.O) GC TO 8Q0 
C 

DO i 1=1, KK 

1 F (P( n .LL.C . ) GO TO 1 
L2=i.-2.*P( I ) /RL+ ( P ( I )/RL/CFPA I **2 
IF (02.LT.0. ) GO TO 1 
t:L = SOKr ( c2 I 
AL=F( 1 » /( 1.-02) 

CF2=( P( 1 )-R2) /i.L/P2 
SF2=-S0i<r ( i,-CF2*CF2 ) 

CF0CK = k 2*(1. + IL*CF2 )-Rl*( 1.+EL*(C12*CF2 + S12*SF2) ) 

IF (ASS (CH..CK ) .GT. 1. ) GO TO 1 
1=2 

F 2= AiNGLr. ( A TAN 2 ( S F 2 , C F 2 ) *R D » 

F 1=F2-ANG12 

FPAL1 = ASIN(l.L*S 1N(F1*0R)/S0RT(1 .+2.*OL*COS( F1*GR) + :L*:L))*RD 
VL1 = S0RT(U*(2./R1-1 ./AD) 

DtL\/=SCRr 1 VC*VG+VL1*VL1-2.*V0*VL1*C0S ( (FPAG-FPAL 1 )*0R ) ) 
STB=VLl/OfcLV*SlN( ( F PAO-FP AL 1 ) *D R ) 

CTh=( V0*V0+DELV*D2L V-VLl*VLl)/2./UcLV/V0 
the f A= AT AN2 (STF ,CTh )*RD 
1 CONTI NU;_- 
C 

IF ( DuLV.LT. .CCl ) GC TC 800 

FL0=F1 

FLC=F2 

K(<=1 

3o0 PtIURN 
lNC 
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SUeROUT CRCSS (XI ,Y 1 , Z1 , X2 . Y2 , Z2 » PX , PY f P 2 , PRODUCT ) 
C CALCULATE VuCrCR CPCSS PRODUCT 

PX=Y1*Z2-Z1*Y2 
PY=Z1*X2-X1>KZ2 
PZ = X1*Y2-Yl=f X2 

PRCDUC f = SORT(PX*PX + PYX=PY + PZ*PZ » 

Ri TURN 
t.ND 


SUeROUriN^ CUEIC (A.B.CtCt X1,X2»X3,KK) 

This SU3R0UTIN£ solves the EOUATION AX##3 +BX**2 +CX +D = 0 FOR 
The R..AL kOCTS 

A.SfC.O - CL^FFICItkT CF Tri*. DIFFERcNT POWERS OF X 
X1,X2,X3 - REAL RJCfS OF THt eOUATION 
KK - NU-10ER OF REAL ROOTS 

CER f( X)=SIGN (ABSI X ) **, 333 333333, X) 

KK = 0 

Pl=3. 1A15927 

I F( A.LF ..1:_-3C . Ai\iC. A. GT.-. IE-3D ) GO TO 4 
P=E/A 
C = C/A 
R=0/A 

SA= (3,*0-P=i'*2)/3. 

S3= ( 2. =1*?* *3- 9. *P#0+ 2 7, *R) 7 27 . 

CcL=( n. *D*»3-C*«2«PY*2-18, *C*P*R+27.#R**2+4,*P**3*R) /1C8. 
IF(Ci.L.LT,.U-3o.ANO.DrL.GT.-.le-30) GO TO 3 
Ih(OeL) 1,3,2 

1 KK = 3 

CPni=-S3/2./S0KT(SA**3/(-27. ) > 

IF (ABS(CPhI ).GT.l, )GO TO 10 
SPhi=SCRT ( l.-CPHl**2) 

P,HI=A FAN2( SPt.I ,CPHl ) 

GO TO 11 

1C SPhI=SOR r { ( 27.*Ll L ) /SA**3 ) 

C SINCe FOK SMALL ANGLES SPhI=PH1 

BL TA = 5PHl 

IF (-SH.GT.C . )Fhl=BclA 

1 F(-SB.L r. : . )P,-iI=3. 1415926535B9793-BCT A 
11 eO=2.*SUKT (-SA/3. ) 

XI =cC*COS(Ph 1/3. )-P/3. 

X2=t:U*C0S( PhI/3.+2.=»PI/3. J-P/3. 

X3 = CU*CfJS( Ph I/3,+4,=»Pl/3. )-P/3, 

GO fO 7 

2 KK =1 

Xl = Cf3R n- SB/2. + SWRT (DLL » > +CBRr ( -SB/2.-S0RT ( DEL »)-P/3. 

GO TO 7 

3 KK = 3 

X 1 = 2 .*CBR r(-SB/2. )-P/3. 
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X2=CBRT(SB/2. )-F/3. 

X3 = X2 
GO TO 7 

4 CCMINUE- 

DI S=C*^2-4,^6=«^D 
IF(DIS) 7»5,5 

5 X1=(-C+S0RT (DIS ) )/2 ./B 
X2=(-C-S0RT (DIS ) ) /2 ./B 
KK=2 

7 CONfINUE 
RfcTUKK 
tND 


SUBROUTINL CtTER(A.e> 

OOtBLi: PR^CISICN APrBP 
DIHENSICN B(3f3) ,BP(3,3) 

DO 10 1=1 f 3 
DO 10 J=l»3 
BP(I,J)=O.DO 
BP(I,J)=DBLt (B( I,J) ) 

10 CCMlNUb 

AP = BP(lti)=i'BP(2,2)*BP(3,3)-8P<3f 1) «BP < 2 » 2 ) «BP ( 1 r 3 * + 
1BP(1,2)^BP(2,3)^BP(3,1)“BP( i , 2 BP ( 2 , 1 ) ^6 P ( 3 , 3 I +8 P (1 , 3 ) *BP ( 2 , 1 ) 
2*BP {3t2)-BP ( 1 , (2f3)=<=BP(3,2) 

A=SNGL( AP) 

RETURN 

END 


SUBROUriNt DOT ( X 1 » Y 1 , 7 1 ♦ X2 , Y 2 , Z 2 , ANGL b ) 

This SUBROUTINE CCN^PUTES THE ANGLE BETWEEN TWO VECTORS 

XltYltZl - CCMPONirNTS GF THE VECTOR R1 
X2,Y2tZ2 - COMPONENTS OF THE VECTOR R2 
ANGLE “ ANGLE BETWEEN VECTORS R1 AND R2 

RD=57,29577S513C823 
R1=S0RT( X1#X1+Y1«Y1+Z1^Z1 I 
R2 = S0RT (X2«X2 + Y2*Y2 +Z2«Z2 ) 

ANGLE = ACOS( ( XI X2 + Y 1 * Y2 +Z 1 « Z2 ) / R i/R 2 ) »RD 

R£ TURN 

END 
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SUBKOUriNc tt4RTH( JD, XhEtYHLf ZHL »DXHL»OYHc»DZHH) 

C 

C This SUSRuUTINE CGVPUTfcS THU hL L I OC cNT RI C PCSITION AND VELOCITY OF 

C The EARTH IN MLAN EQUINOX AND ECLIPTIC OF DATE COORDINATE SYSTEM. 

C This ROUTINE CALLS SUBROUTINES TINVS AND CONCAR, 

C 

C JD - JULIAN DATE 

C Xhc.Yht.Zht - PCSITION OF EARTH 

C CXt-cCYHE.CZhE - VELOCITY OF EARTH 

C 

FOAL JD 

ANGLc (X)=AMCD (X,36C , ) +180 .-S I GN ( 180 . , X ) 

CR=,017A5)329251S9A3 
R0=S7.29b779513C823 
LSLN=1. 3271^ llt+1 1 
AU=1493968A3, 

C 

C=JC-2A15L20, 

C0=C/1C 000. 

Tt: = C/ 3o525. 

C 

A3 = l.C0::O023“!‘AU 

i.h=.;).61o75lC4-0.w J0 04180*T3-O.U0000012 6*Tt:*=i<2 
XI £=0.0 

Ri = i01. 22 J8 3 3+0. ODC C4 7o68 #D+0. 00003 39 >i‘CD**2 +0 . OOGOOOO 7*CD*Xt3 

Oi.=0.c 

XMfc = ANGLt ( 353.47 5 64 5 + 0.98 360026 7*D-O.Oa00112*CO**2-O.CC G0CQC7*C0** 
13) 

C 

CALL TINVSI XML*CR,fc.“,r:Cc,FE) 

CALL CCNCARIAfc ,t2.XIE,WL.0E,FE*RJ,XHE.YHE.ZhL,DXHEiDYH£,DZHE.USUN) 
C 

Ft TURN 
LND 


SUfiROUTlNc LMARSIJO , XHM . Yh M , Z hM , DX rtM , 0 YnM , D ZHM ) 

THIS SUBROUTINc COMPUTlS THc MEAN HELIOCENTRIC POSITION AND 
VELOCITY OF MARS IN THt. MEAN EARTH L'OUINOX AND ECLIPTIC OF DATE 
COORDINATE SYSTEM. THIS ROUTINE CALLS SUBROUTINES TINVS AND CONCAR 

JD - JULIAN DATl 

XhM.YhM.ZnM - PCSITION OF MARS 

CXhM.DYHM.DZHM - VELOCITY OF MARS 


RtAL JD 

ANGLE < X)=AMCO( X , 360 . I +180 .-SIGN ( 180 . t X ) 

0R = .017<+53292519943 

PO=57. 2957795130823 

USUN=I. 327154451+11 

AU=1495988a5. 
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D=JD-2415C20. 

CD=D/10j00. 

Tt=C/36525, 

C 

AM = 1.523(j915*AU 

c«=C.C933129C+O.COCC92C64»<T£-O.OaJOu0077^TE-**2 
Xiy = 1.8503 34-O.COU5 75*Tt+U.000012#TE*’f=2 

CM=48.786442+C.77C991*Tfc-O.Ou30015*rL**2-0,COC00576*T2**3 
wM = 334.218 2C3+l.e4C759*Tt: + O.0OO130*T2#*2-O.COOJQ129»T = **3-OM 
XMM=ANGL=(319.5 2 942 5+3.524020 7666*D+0.OJOO13 35 3*CO**2+0 .0C0C03025* 
1CD**3 ) 

C 

CALL riNVS( XMM*DR,£y ,tCy»FM) 

CALL CCNCAR( AM,EM,XIM,WM,OM,FM*RDtXhM, YHM,ZhM,OXhM,OYHM,DZHM,USUN) 
C 

Rt TURN 
END 


SUBROUTINE LULER < X , Y , Z f XP ,YP , ZP . PH 1 , PS I ,THET A , DPril » DPS I ,DTHETA,WXP 
1,WYP,WZP,J,K) 

XPHI = PHI «,017453292 5 
XPSI=PSI*. 0174532925 
xrh = Th2rA*, Cl 7 45329 25 
lF(J)10tl2,ll 

10 X= (CUSI XPSn*CCS( XFHl )-C0S (XThM^SINIXPHI >*SIN( XPSI ) )*XP + (-SIN( XPSI 
1 )*CUS( XPHl )-CCS (XTH >*SIN( XPHI »»C0S(XPS1 ) »*YP + ( SIN( XTH > * S IN ( XP:U ) )* 
2 IP 

Y=(COS(XPSn#SlN(XFHl » +COS ( X TH ) *C0S ( XPH I )#SIN( XPSI » > *XP+ (-S 1 N I XPS 1 
U4SIN(XPHi)+CCS( XTh l*COS( XPMl I’XCUSI XPSI ) ) <-YP + ( -SIN ( XTH I *COS ( XPHl ) ) 
2«ZP 

Z= ( SIM XTh)»SIN( XPS I ) »*XP + ( SIN( XTH I *CUS( XPSn ) *YP+ (C0S( XTH » >*ZP 
GO TO 12 

1 1 XP=(CUS( XPSn *C0S( XPnl »-CUS( XTHM=SiN<XPhI |4SIN(XPSI >) ■S‘X+ (CDS ( XPSI ) 
1*S1M XPHl )+CUS( XTH) ^CCSIXPHl )*S I N ( X PS 1 )) *Y+ { S I N ( X TH > *S I N ( XPS I )) *Z 

YP=(-S1N( XPSI )*CLS( XPHI )-COS(XTh)*SIN( XPHI )*CUS( XPSI))*X+(-SIN(XPS 
1 1 ) «SIN(XPhI )+C0S(XTH)*CCS ( XPHI )*COS( XPSI ) )*Y+IC0S( XPSI ) *SIN( XTH) )* 
21 

ZP= (SIN (XTH IX-SIM XPHI ) )*X+(-SIN (XTh)*COS( XPHI ) )*Y + C0S( XTH) *Z 

12 1F(K)13,15.14 

13 DPhI=(WXP*S INI XPSI ) +WYP«CUS( XPSI ) )/SlN( XTH) 

CPSI=WZP-(CCS(XTh)* (WXP*SIN(XPSI )+WYP*COS( XPSI ) ) )/SIN( XTm) 

uri-ir-TA=v<xp^=cos(xpsi )-wyp*sin( xpsi ) 

GO TO 15 

14 WXP=DPHI*SIM XTH )*S IN (XPSI )*DTHE TA4C0S ( XPS I ) 

WYP = DPHI*SIN(XTF )*CCS (XPSI )-DTHL TAX'S in ( XPSI ) 

KZP=DPSI+OPhIX‘COS(XT)n) 

15 RETURN 
END 
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SUBRQUTINl JULCAK XfWDI iFUI f INO) 

This SUBRUUTINa CCMVLRTS A GIVEN JULIAN DATE OR THE NUMBER OF 
WHCLC AND FRACTIONAL CAYS SINCE JANUARY 1, 1950, 0 HRS*, TO THE 

CUKRESPONDING CALENCAR DATE. 

WDI - integral part of JULIAN DATE OR WHOLE NUMBER OF DAYS SINCE 
JANUARY 1, 195Cf G HRS* 

FDI - FRACTIONAL PART OF JULIAN DATE OR FRACTIONAL NUMBER OF DAYS 
SINCE JANUARY 1, 195Ut 0 HRS. 

INO - CONTROL INTEGcR. 0 IMPLIcS JULIAN DATt,l INPLIES DAYS 
X(l-b) - CALENDAR CATE ( YEAR , MONTH , DAY , HOUR , MI NUTE , SECOND ) 

DIMENSION X (6) ,A(12 ) ,W(12 ) 

wG = WD 1 

FC=FD1 

lF(IND)ltl,5 

1 IF (FD-.5) 2*2,3 

2 FD=FD+.5 
WC = v^D-l. 

GO TO A 

3 FC=FD-.5 

A WC=vmD--2a33282. 

5 WC=WD+1. 

CY=365. 

Z = 2. 

N='J 
0 = A. 

6 WC=WD"DY 
IF(WD) 10,10.7 

7 N=K+1 
Z=Z+l. 

CK-C-Z 

IF(CK)9»9.B 

8 CY=3o5. 

FC=28. 

GO TO 6 

9 0Y=366. 

Q=C+4. 

FC=29. 

GO TO 6 

10 WO=i^D + CY 

DU li 1=1.12 

11 A( i ) = 0. 

Cl = 31 • 

C2=3G. 

DO 13 1=1,12 
A( I ) = 1 . 

CA=FC^A( 2)+C1^(A(1)+A(3)+A(5)+A(7)+A(S)+A(1C)+A(12))+C2«(A(4)+ 

1 A(6)+a(9)+A( 11) ) 

W ( I ) = wC«CA 
1F(W( I ) ) 12,12, 13 

12 1F( 1-1) 15,15,16 
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15 MGN=1 

GO TO 14 

16 MGN = I-1 
fc^D = W(MCN.) 

^CN=MCN+1 
GO rO 14 

13 CCMINU^ 

1^ N=N+5U 
X( 1 )=N 
X( 2)=M0N 
X(3)=WD 
FH=F0«24. 

N=FH 
X(4) = IM 

FM=(FH-X(4) )*60. 

N = FM 
X( 5)=N 

X(6)=(FM-*X( 5) )«60. 

Ft TURN 

fcNO 


SUBROUTINL LARBRT<Xl,Yi,Zl,X2,Y2,Z2tTIMci2»AtE,Xl , W , 0, TAl , TA2 t U tKK 
1) 

IF X2,Y2,AND Z2 AR£ ZtROt ThciNI XI IS CONSlOtRLO Trie KAOIAL 
DISTANCt TO POINT 1, Y1 THb UISTANCl TO POINT 2, AND Z1 THL 
ANGLt FRGRi POINT 1 TO POINT 2. MOTION IS ALWAYS C0NS1ULR2D 
FROM POINT 1 TO POINT 2. 

RtAL M12fN 

ATANH( X)=.5^AL0G( ( 1 . + X ) / ( 1 •- X M 

ANGLt ( X)=AMCD( X»360 . ) +180 .-S I GN ( 180* »X ) 

CATA CR • RD • PI/. 01 7^ 5 329 25 IS 94, 57. 29 57795 ♦ 3. 141592653 5/ 

M=0 

KK=l 

TA2=0. 

KFY = 1 
IGR8IT=1 

IF ( TIKl12.L£ .0. ) GO TO 8DU 

IF<(A6S(X2).LT..l).ANC.(ABS(Y2).LT..l) .AND. ( AB S ( Z 2 ) . L T . . 1 ) ) GO TO 1 
K1 = SQKT (X1^X1 + Y1*^Y1+Z1^ZU 
R2 = SQK T( X2«X2 + Y2*Y2 + Z2*Z2 ) 

CPSI= (X1<^X2 + Y1«Y2 + Z 1*Z2)/K1/R2 

SPSI = (X1«Y2-X2^Y1 )/ABS(Xl=«^Y2-X2*Yi)=<'S0Kr (l.-CPSI^i'CPSI ) 

PSI=ANGLL( ATAN2 (SPS 1, CPSI i *RD) 

GO TO 2 

1 Rl^Xl 
R2 = Y1 

PSI = ANGL£( Z1 ) 

XI=C. 

0=0. 

Vi = 0. 

2 C=SCRT(R1*R1 + R2X^R2-2.*R1*R2»C0S ( PSI*DR ) ) 
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lF(PSI,Lr.C .01 » GO TC 800 
AM=(Ri+R2+C »/4, 

S=2,#ANi 

rP=S0Rr(2./U*(S<'*1.5-(S-C>**l.5»/3. 
TTP=S0RT(2./U)*(S=!'*1.5+(S-C»*=f=1.5»/3. 
iP( (PSI.L£..180.),ANO.(TIMt:12.Lr.TP) ) 1 0KB 11=2 
IF (( PSI. Gb. 130. ),ANC. (I IME12 .lt, TTPM 10R8I T=2 

3 CTA2=CGS( TA2*0R) 

CrAl=C0S( (TA2-PS1 >*CR> 

0=R2*CTA2-R1*C TAl 
lF(AtSS(0>.GT.l.) GC TO 5 

4 IF(KEY.GT,1)GC TO 2 li 
TA2=TA2+b, 

GO fO 3 
a F=(Rl-K2l/0 

1F(:;,LT,0,» GC TC 4 
A=H2*( l.+„*CTA2 )/< 1 .-b*L) 

GO TO (6,7) .ICRBIT 

6 lF(tl,GT.l..CR.A.LT.O. ) GO TO 4 
TtMP = S(JRT( (1,-L)/(1. + U» 

tCl=ANGLE (2.*ATAM( TcMP*TAN( (TA2-PS1 )*DR/2.) )*RO)*OR 
LC2 = AfMGLb(2.*ATAM( T£MP*TAN( TA2*DR/2, ) ) =!'RD)*DR 
OS L S C = :;C 2— cC 1 

IFCCsiLlC.Lf ,C, ) DLL5C = 2.*P I+DtLtC 
M12 = U£;LtC-S*(SIM(bC2)-SlM(ECl ) ) 

Gu ro 3 

7 lF(t,.LT.l..CR. A.GT. J. ) GO TO 4 
TcMP=SORT( (t-1. )/(fc+l. ) ) 

cCl = 2.*ArAN)-( TuRP*TAN( ( TA2-PS1) #0R/2, ) ) 
uC2=2.*ATANT ( T£MP*TAN( (TA2 *l)R/2. ) ) ) 

M12=.'.>l‘( SiN) (SC2)-SlNh(FCi) )-iC2 + cCl 

8 N=S0kr(U/ABS(A)**3) 

F=TIMtl2-M12/N 

GO TO (9,10,11) ,KtY 

9 K£Y=2 
TALaST=TA2 
TA2=TA2+1. 

GO TO 13 

10 K£Y=3 

GO TO 12 

11 M=M+1 

IF(M.GT.6C)G0 TC 8C J 
IF(AHS(F),LL.ABS(FLAST) ) GO TO 12 
25 DFCTA2=0FCTA2*2. 

M = R+1 

IF(M.GT.6i, )G0 TC 80 0 
TA2= TALAS r-F LA ST/UFDTA2 
GO TO 3 

12 bRPCR=F/TIMtI2 
lF(ABS(i-RRUP),LT. ,00001) GO TO 14 
UFCTA2 =( F-F last )/( TA2-T ALAS T) 

Talas t=ta 2 
TA2=TA2~F/GFCrA2 
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13 FLAST=F 
GO TO 3 

14 TA1=TA2-PSI 

IF((ABS(X2)*LT#*1)#AND# ( AB S ( Y2 I . LT. . 1) • A, ( ABS ( Z2 ) . LT • • 1 ) )GO TO 900 

D1=Y1*Z2-Z1*Y2 

02=Zi*X2-Xl*Z2 

D3=X1*Y2-Y1*X2 

HH=S0RT(Dl=«'Dl + 02*02+D3«D3) 

IF(D3.GT.O.) GO TO 15 

D1=-D1 

02=-D2 

C3=-03 

15 C0SXI=D3/HH 

XI=ATAN2( SORT (l.-COSXI*COSXI ) tCOSXI )«R0 
S0=D1/ (Hii^SINC XI^DR ) ) 

CO=-D2/(HH*SIN(XI’4‘DR) ) 

IF ( SO, £0.0. .AND .CO. £0.0. ) C0=1. 

0=ANGLr C ATAN2(SCtCC)««D) 

W=ANGLfc( ATAN2 ( (-X1«S0 + Y1*C0 ) *COSXI +Z1*SIN ( X I«DR ) t Xl«CO + Y 1*S0 ) »RD-T 
lAl ) 

GO TO 900 
8C0 KK = 0 
900 RETURN 
END 


SUBROUTINii LATLNGC X » Y , Z , XL AT t XLUNG ) 

THIS SUBROUTINE CC^PUTLS THE LATITUDE AND LONGITUOL OF A GlVLN 
POSITION VECTOR 

X*Y*Z -- CCNFCNZNTS OF ThC POSITION VtCTOR 
XLATtXLUNG - LATITUDE AND LONGITUDl 

ARCOS(X)=ACCS( X) 

ARSIN (X) =AS IN( X ) 

RD=57. 2957795 

P=SORT (X«>J'2 + Y*>^2 + Z*«2) 

XLGNG = ATAN2 ( Y,X J<^RC 
XLAT=ARSIN( Z/R ) «RD 
RETURN 
END 
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SOBi^OUTINc. OCCULU A,c»XI f WtOfUf RS.ex ,uY,£Z,TOCC,Tl»ALTl»Fl »DECltR 
lAl t r2.ALT2,F2.0£C2,RA2tKK> 

THIS SUBROUTINE COMFUTtS THE ENTRANCE AND EXIT TRUE ANOMALIES OF 
GCCULTATION. THIS ROUTINE CALLS SUBROUTINES RXYZPQW, QARTIC. 

CROSS, OOr, TCCMC, RPCWXYZ, AND LATLNG. 

A,L,XI - semimajor AXIS, ECCENTRICITY, INCLINATION 
W,0 - ARGUMENT OF PERIAPSIS, LONGITUDE OF ASCENDING NODE 
U,RS - GRAVITATIONAL CONSTANT AND RADIUS OF THE PLANET 
tX,cY,EZ - CGf'PCNENTS OF UNIT VECTOR TOWARD THE BODY OCCULTED 
TOCC - LENGTH OF TIME IN SHADOW 

Tl,ALTl,Fl,DcCl,RAI - CONDITIONS AT ENTRY INTO THE SHADOW, TIME 
FROM PERIAPSIS, ALTITUDE, TRUu ANOMALY, DECLINATION, AND 
RIGHT ASCENSION 

T2,ALT2,F2,CEC2,RA2 - CCNDITIUNS AT LXIT FROM THE SHADOW 

KK - CCNTROL INTcGtR. 0 IMPLIES NO OCCULTATION, 1 IMPLIES OCCULT. 

DIMENSICN RP0W(3,3) ,CF (4) , RXY2 ( 3 , 3 ) 

ANGLt (X»=AMCD(X, 360 .) +180. -SIGN ( 180. ,X> 

CR=. 017453292519943 
PO=57.2S57795130823 
Fl = 2( v^O. 

F2=20C0. 

KK = 0 

P=A*( l-i*£r) 

CALL RXYZPOW(rX,tY,EZ,XI , W , 0 , RPOW , BE T A , XXI ,ZBODY» 

Cl = {RS/P>**4*E**4-2.*(RS/P)**2*(XXI«*2-8fcTA**2»X'E**2 + (B£TA**2 + XXI* 
1*2 ) **2 

C2 = 4.*(RS/P l**4*L** 3-4.*(RS/P) **2*( XX1**2-BETA**2»*E 
C3 = 6.*(RS/PJ**4*L**2-2.*( P S/P I * *2* ( XX I **2-BETA**2 I -2 . * ( RS/ P » **2*( 1 
1,-XXI**2 )*E**2+2.*( XXI**2-BETA**2»*( 1 . -XXI **2 > -4. *BETA**2* XXI **2 
C4=4.* (RS/P )**4*E-4,*(RS/P»**2*( l.-XXI**2»*E 
C5=(RS/P)**4-2.*(RS/P )**2*( 1 XX I **2 ) + ( 1 . -X XI **2 » **2 
CALL yARTIC(Cl,C2,C3,C4,C5,CF( U , CF ( 2 ) , CF ( 3 » ,CF ( 4 > , J J ) 

IFUJ.tO.OI GO TC 8C0 

DU 3 1=1, JJ 

IF (A8S(CF( I n .LT.l.OC Jl.AND.ABS (CF( I n .GT, 1. ) CF { H =0. 999999 
IF (ABS(CF( n ).GT.l. ) GO TO 3 
SF=SOHn i.-CF( I )**2 I 
R = P/(1.+E*CF(I) ) 

CALL CROSS (R*CF ( I ) , K* SF ,0 , , BE TA , XX I , Z8UDY , PX , PY , PZ , PRODUCT ) 

CALL DOriCFd) , SF, D., BETA, XXI, Z BODY, ANGI 
IF(ABS{PR0DUCT-PS).LT.,01.AND,ANG.GT.90.) GO TO 1 
SF=-SF 

CALL CRIJSS(R*CF(I) , P* SF ,0. , BETA , XXI , ZBODY, P X , PY, PZ, PRODUCT ) 

CALL CCTICFU ) , SF , C . , BETA , XX I , ZBODY , AN G ) 
IF(ABS(PRODUCT-RS>.LT,.01.AND.ANG.GT.90.) GO TO 1 
GO TO 3 

1 IF(F1.LT,10C0. >G0 TO 2 

F1 = ANGLc(ATAN2(SF,CF( I) )*RD) 

GC ro 3 

2 F2 = ANGLh(ArAN2(SF,CF(in*R0J 
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GO TQ 4 

3 CONTINUE 
C 

4 IF(F2.GT.1CC0. )GO TQ 300 
CF1 = CUS(F1’4'CR) 

SF1=S1N(F1^DR) 

DSDF=2.^RS’«'RS*( l.+fc^CFl)^(-t*SFl ) + 2, ( 8E TA^CF 1 + XX I « SFl )=4^(-BETA 

l«SFi+XXI«CFl) 

1F(CSDF.GT.C. )G0 TO 5 
C 

TtMP=Fl 
F1 = F2 
F2=FEMP 
C 

5 CALL TCONIC(U»c,A,Fl,Tl) 

CALL TCQNIC(U,e,A»F2,T2) 

T0CC=T2-ri 

IF (TOCC.LT.O. ) TUCC = 6.26 3ie5 3^S0Rn A«*3/U)+T0CC 
Tl=Tl/60. 

T2=T2/60, 

T0CC=rUCC/6C, 

ALTl = P/< 1,+E*C0S(F1=4'CR) )-RS 
ALT2=P/(1.+L«CCS(F2^DR) )-RS 

CALL RPOWXYZ ( C0S(F1«DR),S1N(F1^DR>,0.*X1,W,0,RXYZ,RX,RY,RZ) 

CALL LArLNG(RX»RY,RZ»OtClfRAli 

CALL RP0WXYZ(CGS(F2^DR) rSIN(F2*DK) ,0. ,XI ,W,0,RXYZ,RX,RY,RZJ 
CALL LArLNG(RX»PY»R 2» CLC2*RA2) 

KK = 1 

GO rO 930 
C 

SCO CCNTINUt 
T0CC=0, 

T1=0. 

ALT1=0. 

F1 = 0. 

ObCl-0. 

RA1=0. 

T2 = 0. 

AL I2=J. 

F2=0. 

CfcC2=0« 

RA2=0. 

900 CONTINUE 
RETURN 
END 
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SUHROUriNE PARINCXP ,YPtX, Y, in 
DINENSIGN X(3) fY(3) fOET(3»3l 
DO 1C 1=1.3 
OfcTU.l )=X( n*«2 
DLT( I.2)=X(n 
DcT( I.3> = i,C 
10 CCMINUc 

CALL DLTLR(D,CLTI 
1F(A6S(D>-1.0E-16) 100,100,20 
20 DO 30 1=1,3 
DET( l,l)=Y(n 
30 CCMINUH 

CALL DfcFLRt A,DET) 

A = A/D 

CO 4j 1=1,3 
DcT(I.l)=X(n*=«'2 
CcT ( I.2)=Y(I ) 

40 CONTINUE 

CALL DcTuRiBfOLl) 

B=8/0 

DO 50 1=1,3 
DLT(I.2)=X(I) 

Ct T( I ,3)=Y( I ) 

50 CONTINUE 

CALL CtTLR(C,DLn 
C=C/D 

lF(in2C0,15C,2Cj 
150 XP=-B/(2. Q*A) 

YP = L-[3^«2/(4,C*A) 

GO TO 3lO 

ZOO YP=(A-^XP + B)^XP + C 
Gu ra 300 
ICC XP=C.G 
YP=0.0 
3C0 RcTURN 
tNC 


SUBROUTIN:- PRLCLS( JJI ,XE1,YE1, ZL1,J02,X£2,Y£2,ZE2> 

THIS SUBROUTINE TRANSFORMS GLUCfcNTRIC EARTH EQUATORIAL CGOROINATtS 
FPCv lPOCH JDl TL EPOCH JD2. THIS ROUTINE CALLS SUBROUTINE EULER. 

JD1.JD2 - JULIAN DATES OF INITIAL AND FINAL EPOCH 
XLl.YLlfZEl - COMPONENTS OF VECTOR IN JOl COORDINATE SYSTEM 
XL2 .Yl2,ZE2 - COMPONENTS OF VECTOR IN JD2 COORDINATE SYSTEM 

REAL JC1,JD2 

T=ABS( ( JD2-JD1 )/36524.219679) 

TC= (JD2-241 5020. ) 736524. 2 19879 

ZETA0=(0.64CC6944^C .38777 778E-3* TO) ’5'T + O.83888839e-4*T«»2+0.5E-5«T^ 
1«3 
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C^ETA0=ZETAC+C,21972222t-3«r«*2 

THETA0=<0. t>568 5611-3. 2369 AA4fc- 3«T0 ) *T-0. 1 1 8 33333E-3*T«* 2-0 • 1 1 66666 
17c-4«T>i^^3 
C 

IF( JD2-JO1.GT.0.) GO TO 1 
TEMP=ZETAO 
ZETAC=-CZCTA0 
CZETA0=-rtMP 
ThE7Au=-THtTA0 
C 

1 CALL £ULcR(XElfYtl»Ztl,Xt:2»YE2»Z32,90,-ZtTA0,-(90,+CZtTA0) ,TH3TA0f 
IDPhI »UPSl,OPSI»V^XP,^YP,^<ZP#l,0) 

C 

RETURN 

END 


SUBROUTINE CACRAT<A»BfC.XltX2fKK) 

SOLVES EOUATIGN A^X«*2 + B*X+C=0 
KK = NUMBER OF REAL ROOTS 

KK = C 
DI S = 

IF (DIS.LT.O. ) GO TO 800 
Xl=(-B+SgRT (DIS) 

X2=(-B-S0KT(DIS) >/2./A 
KK = 2 

BOO RETURN 
END 


SUBROUTINE C AR T I C ( A t B ♦ C , 0 , E t X 1 , X2 f X3»X4»KK> 

THIS SUBROUTINE SOLVES THE EQUATION AX**4 +BX^*3 +CX««2 +0X +E ^ 0 
FOR THE REAL RUCTS* THIS ROUTINE CALLS SUBROUTINcS QADRAT AND CONIC 

A.BfCtDtb - CCEFFICILNTS OF THE DIFFERENT PCWLRS OF X 
Xl,X2fX3,X4 - REAL RLGTS OF THE EQUATION 
KK - NUMBi-R CF REAL ROOTS 

KK = C 
BP=B/A 
CP=C/A 
CP=C/A 
LP=6/A 
C 

H=-8P/4. 

H 2 =H **2 
h3=H2*H 
F»4 = Fj3«H 

P = e.*H2 + 3.«BP’«'H + CP 

G=4.^h3+3.«BP^H2<'2,^CP«H+0P 

R=H4+BP*H3+CP«H2+DP#h+£P 
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CALL CUBIC ( 1. t2.’{'P,P*P-A,^R,-Q=<'0,Tl,T2,T3»NR00n 

GO TO ( 1,2.3 ) .NROOr 

1 RP=T1 
GO TO 4 

2 RP=AMAxi(TlfT2l 
GG TO 4 

3 RP=AMAXi(Tl,T2,T3) 

4 CCMINUL 
SORP=SGRT(RP) 

XI=(P+RP-0/SQkP)/2, 

RcTA=(P+RP+C/SGRP)/2, 

CALL QADRAT (1. fSCRP ,XI,Y1, Y2, IROQT) 
call CADRATd* ,-SCHP, 3LTAfY3,Y4, JROOT ) 

1F{ iRQCT + JRCiOT.tO.O ) GO TO 800 

IF(lKCCr+JRCCr.£Q,4) GO TO 6 

IF( IRUGT.tQ.C ) GC TC 5 

X1=Y1+H 

X2=Y2+h 

KK = 2 

GO r 0 s 00 

5 Xl=Y3+ri 
X2=Y4+r 
KK = 2 

GC ro 800 

6 Xl=Yl+h 
X2=Y2+n 
X3=Y3+h 
X4=Y4+ii 
KK = 4 

8v.O CCKTIMIr 
Pc rUHM 
LNC 


SUePOUriNt RcCCCC JO ,X£C.Y cC.Z£C,X£OtYPO,ZcO) 

This SUBRuU^I^5 P0T4T£S A VcCfOR FROM G£OCtNrRIC, tCLIPTIC. TO 
THE GcOCEMrPIC. cARTh tQUATGRIAL COORDINATE SYSTEM 

JD - JULIAN DATE 

XECtYLurZuC - CCMPCNLNTS OF THE VECTOR IN THE GEOCENTRIC, ECLIPTIC 
COORDINATE SYSTEM 

XfcC,YcQ,ZEO - COMPONENTS OF THE VECTOR IN THE GEOCENTRIC, EARTH 
EfcUATCRIAL, COORDINATE SYSTEM 


RlAL jc 

OR = .J17A53292!3l9943 
Tt=( JC-241bC20. )/36525. 

XIE=23.A52294-O.G13C125*TE-O.OC000164*Te**2+O.OOOOQ0503*TE**3 
C=CGS(Xlfc*OR ) 

S=SIN(XIE>!'DR> 
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XEC=XEC 

V£C=YEC*C-ZfcC*S 

ZEC=YEC*S+Z£C*C 

RETURN 

END 


SUBROUTINE RE0ME0( JC» XEO f YcQ , ZEO t XMEQ , YMEO . ZMEQ , DcCMcQ , R AMEO ) 

C 

C THIS SUBROUTINE ROTATES A VECTOR FROM THE MEAN EARTH EQUATOR- 

C EQUINOX TO THE MEAN MAPS EQUATOR-EQUINOX COORDINATE SYSTEM. THIS 

C ROUTINE CALLS SUBROUTINES REQPcQ AND LATLNG. 

C 

C JD - JULIAN DATE AT TIME OF INTEREST 

C XEG.YEO.ZcO - VECTOR IN THE EARTH EQUATORIAL SYSTEM 

C XMECtYMEQ.ZMEC - VECTOR IN THE MARS EQUATORIAL SYSTEM 

C DECMcC.RAMEQ - CECLINATION AND RIGHT ASCENSION OF THE VECTOR IN 

C THE MARS EQUATORIAL SYSTEM 

C 

REAL JD 
C 

TL=(JD-2A15020. )/36525. 

TAU=AMQO{ TE*1J0. .1. ) 

TP=TE“«'100.-TAU-50. 

C 

ALPHAO = 316.5 5+4 5. OC 6751 + . 006 7 Sl-S'TP- . 00 101 3#TAU 
GAMMA0 = 52. 8 5+A5*.CiL 343 + .003A8«TP-.000631#TAU 
OMEGA=48.78£44167+0.77C99167*TE-0. 13888889L-5*TF**2 
XI = 1.8 503333 33-C.67 5E-3=f‘Te+0,12611111E-4#TE*«2 
C 

CALL REQPEOIJCi ALPH AO , GAMMAO , OMEGA , XI , XEO .YE C f ZEO . XMEQ , YMEO t ZMEQ) 
CALL LATLNGI XMEC, YMEQ, ZMEQ .DECMEQ.RAMLO) 

C 

RE TURN 
END 
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SUBRUUriNc PlCPEO( JCt/lLPHAC.GAMMAOf UMLGAtXI , XEQ , YEQ, ZEQ » XPEO f Y PEO t 
IZPcCI 
C 

C THIS SUBROUTINE ROTATES A VECTOR FROM MEAN EARTH EQUATOR-ECUINOX 

C TO PLANET ECUATCR-ECU INGX COORDINATE SYSTEM, THIS ROUTINE CALLS 

C SUBROUTINE LULER. 

C 

C JU - JULIAN DATE AT TIME OF INTEREST 

C ALPHAG.GAMMAO - RIGHT ASCENSION AND DECLINATION OF THE PLANETS 

C AXIS CF RCTATICN EXPRESSED IN THE EARTH EQUATORIAL 

C COORDINATE SYSTEM 

C OMEGA, XI - LONGITUDE OF THE ASCENDING NODc AND INCLINATION OF THE 

C PLANETS ORBITAL PLANE REFERENCED TO THE ECLIPTIC AND 

C VERNAL kiCUINOX 

C XEC .YiiQ, Z'iO - CCMPCNENTS OF THE VECTOR IN THE EARTH EQUATORIAL 

C COORDINATE SYSTEM 

C XPtC, YPEO,ZPEO - CGMPCNtNTS OF THE VECTOR IN THE PLANET EQUATORIAL 

C COORDINATE SYSTEM 

RcAL JO 

DIMcNSICN RFQW(3,3) 

CR=0. 01745329251994 3 
RD=57.295779513L£23 
TL= ( JC-2415C20. 1/36525, 

E; = 23.4522 9444-C ,13012 5E-1*TE-C, 16388889E-5 X‘Tt*#2+0 , 5D277778E-6*TE * 
1*3 
C 

CL=COS( L*OR ) 

SE = SlN(i.*DR ) 

CAL=COS( ALPHAC*CR) 

SAL=SIN( ALPFAO+DR) 

CGM=C0S(GAMRAC*CRI 

SGM=SIN(GAMRAC*DR) 

CCM=CCS(CMcGA*OR) 

SOR=SIN(OMLGA*OR) 

C 

CZP=CH*SOM*CAL-CCM*SAL 
SZP=SCRT a.-CZP*CZP > 

ZP = ATAN2.(SZP,CZF)*BU 
SXP=Sl.-*CAL/SZP 

CXP = (-Cc*a>*CAL-SCM*SAL)/SZP 

XP=ATAN2(SXP,CXP)*RC 

SYP = Si-*SCR/SZP 

CYP= (Cl*SOM*SAL+CCM*CAL )/SZP 
VP = ATAN2(SYP,CYF)*Rf' 

CI=COS( (XP-XI )*DR)*S1N( (YP-GAMMAQ)*DR»+S1N( ( XP-X I> *0R ) *COS ( ( YP-GAM 
1RAC)*DR»*CZP 

si=soRr(i.-ci*cn 

SWF=SZP*SIN( (XP-X1I*DR)/SI 

CWP=(-CUS( ( XP-XI )*DR) *COS( ( YP-G AMMAO) *0R ) +S IN ( ( XP-XI)*DR»*SIN((YP- 
1GAMMAC)*0R)*CZP)/SI 
WP=ATAN2(SWP,CWF )*RD 

CALL ZULER(Xt0,YEG,Zh0,XPt0»YPE0»ZPEQ,9a,4-ALPhAO,WP + 180. ,90, -GAMMA 
l(.j,L,,0,,O,fC,,0,,U, ,1,0) 

RL TURN 
tND 
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SUB ROUTINE RPQWXYZi VPfVGt VWtXI t W t 0 t RX YZ ♦ VXf V Y , VZ ) 

C 

C THIS SUBROUTINE ROTATES A VECTOR FROM THE PQW TO THE XYZ 

C COORDINATE SYSTEM* 

C 

C VPfVO.VW - CC^FCNE^^S OF THE VcCTUR IN THE PQW SYSTEM 

C XI,W»0 - INCLINAIICN, ARGUMENT OF PEKIAPSISt LONGITUDE OF ASCENDING 

C NOCL 

C RXYZ - ROTATIONAL MATRIX FROM THE POW TO THE XYZ COORDINATE SYSTEM 

C VX.VYfVZ - CGMFCNENTS OF THE VECTOR IN THE XYZ SYSTEM 

C 

CIMcNSIGN RXYZ(3,3) 

CR=.0174532S25199A3 

C 

CW=CGS( W*DR) 

SW = SINC W*OR ) 

CG=CUS(0*0R) 

S0=SIN(0*DR) 

CXI=COS( XI*CR ) 

SXI^SIN( XI^DR) 

C 

RXYZ( 1 »l)=CW*CO-SW*SO=4^CXI 
RXYZ { I ,2 )=-SW*Ca-CW*SO«CXI 
RXYZd »3l = SO*SXI 
RXYZ(2f 1 )=CW«SO+SW*CO^CXI 
RXV Z(2,2)=-SW’!^SO + CW«CQ*CXI 
RXVZ(2,3)=-C0«SXI 
RXYZ(3»1)=SW^SX1 
RXYZ( 3,2)=CW*SXI 
RXYZ<3,3)=CXI 
C 

VX = RXYZ ( !♦ 1 )«VP + RXYZ( 1 , 2 ) ^VQ+'R X Y Z C 1 1 3 ) =«=VW 
VY=RXYZ(2f 1 )«VP+RXYZ(2t2)*VO+RXYZ(2,3I^VW 
VZ = RXYZ( 3, 1 )*VP+RXYZ( 3,2) »VQ + RXYZ< 3,3) *VW 
C 

RETURN 

END 


SURRCUriNL RXYZPQU \*X ,VY, VZ,XI , W ,0 , RPOW , VP , VO , VW ) 

ROfATLS A VECTOR FROM THE XYZ TO THE POW COORDINATE SYSTEM 

DIMENSION RPQW(3,3) 

DR=.01745329251S943 
C 

CW=COS( W^DR ) 

SW=SIN(W^DR) 

C0=C0S(0*DR) 

SO=SINCQ’«'DR) 

CXI=CCS( XI^^'DR) 

SX1=S1N(XI^CR) 

C 
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hP:Jv^{ 1 ,1 )=CW=f'CC~SW*SO*CXI 
PP0^< 1 ,2 ) =Ck>i=SC + SW=)‘CC*CXI 
PPGW( 1 ,3»=SVi*SXl 
PPCW(2,1 )=-SW*CC-Ch*SG*CXI 
RPCV^(2,2 )=-SW=!‘SO + Cw*Ca*CXI 
PP0vJ(2f3) = C^*SXI 
PPCi^(3,l)=SC*SXI 

kPCw ( 3 , 2 ) =-ca*sxi 

RPQW( 3,3 )=CXI 
C 

VP=PP0W(1,1)*VX+PPQW( 1 ,2)-VY+RPQW( 1 , 3 ) *VZ 
VG=PPQW(2,1 )*VX+RPCrt(2,2)*VY+RPQW(2,3) *VZ 
VW=KP0w(3,l)*VX+FPQW( 3 , 2 ) * V Y+R P QW ( 3 , 3 ) * VZ 
C 

PS: TURN 

U.NC 


SUBROUT I Nc TCGiN IC (U , EC , A , TA , T ) 

DATA Ok/. 17A5329251SSR3E-1/ 

T A2=TA*DR 
SLR=A# ( l.~i-C*2C » 

AB = ABS{ A I 
FAC=AB*SORT ( AB/U ) 
cCA=( l.~EC)/(l.+uC) 

AriE = SOKT ( AES ( EG A ) ) 

Th,:=TAN( .3*TA2) 

IF(ABt:-.5t-10m,ll,12 

12 CCNTINUc 
ECA=2.*ATAN(AS2*Thi; ) 

IF(A)14,U,13 

13 T = FAC*(LCA-LC’^SlN(tCA) ) 

GU TO 13 

lA ANG = ABt*T,T2 

ANG = 1.+2.»^ANG/(1.-ANG) 

T=FAC*(cC*TAN(LCA)-ALUG(ANG») 

GO TO 16 

11 FAC = SGk T ( SLR*=(-'3/U) <<2./ ( (l. + tC)**21 
ECl = tCA*TUE *«2 

T = FAC*( Th r + n :*>!»3*( ( 1 .-2 . *LC A ) / 3 .- ( 2 . - 3 . *EC A » *LC 1 / 5 . + ( 3.-4.*ECA)*E 
lCl«^2/7.-(A.-5.*tCA)=i=ECl**3/9,)l 

16 CCNTlNUr. 

RE TURN 
END 
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APPENDIX D — Continued 


SUBROUTINE T I N VS < M ,E ,EC t F ) 

REAL M,MO 

DATA PI/3, 1A1592653589793/ 

ASINh(X).= SIGN« ALOGI ABS(X)+SOKTI X**2+l. ) ) »X) 
IFIt.GL.l. )G0 TO ICJ 
fcC = M 

10 MU=cC-c#SINC=C) 

CM=M-M0 

Dt=CM/(l,-L*CCS(EC) ) 

£C=cC+DE 

lF<ABS(De).GT.l,E-12 )GO TO 10 
HEC= ec/2. 

hF=ATAN(SORT( <l.+U/( 1 ,-H )) *SI N ( HEC ) /COS ( ribC ) ) 
IFIHF.LT.O. )I1F=HF + PI 
F=2.*hF 
GO TO 800 
ICO CONTINUE 

cC=ASlNH(M/E) 

101 KO=£*SINH(EC)-£C 
CM=R-MO 

DE = CM/ (c*COSHIEC)-l , » 
tC=EC+Dt 

IF (ABS(0E>.GT.1,L-12 ) GO TO 101 

F=2.*ATAN( SQRT( (i; + l,C )/<E-l,0» » *TANH( tC/2 .0 ) » 
800 RtfURN 
END 


SUBROUTINE VECTOR ( J D , DECS .RAS . DEC c ,RAE ,DECC ,RAC » SX, SY , S Z ,£X , EY , t Z , 
1CX,CY,CZ,IBCDY> 

This SUBROUTINE CGRPUTtS TF£ POSITION OF THE SUN, EARTH, AND 
CANOPUS IN PLANlT EQUATOR, MEAN PLANET EQUINOX OF DATE, AND WRITES 
DATA. THIS ROUTINE CALLS SUBROUTINES EEARTH, EMARS, EVENUS, PRECES, 
LATLNG, DOT, RECEO, RcGVEO, AND REOMEO, 

JD - JULIAN DATE AT TIME UF INTEREST 

IBODY - CONTROL INTEGER. 2 IMPLIES VENUS, 4 IMPLIES MARS. 

DuCS,RAS - CLCLINATION AND RIGHT ASClVSION OF THE SUN. 

OtCt.KAE - DECLINATION AND RIGHT ASCcNSION OF THE EARTH. 

DECC.RAC - DECLINATION AND RIGHT ASCENSION OF CANOPUS. 

SX,SY,SZ - UNIT VlCTCR FROM THE PLANET TO THE SUN. 

EX,EY,EZ - UNIT VECTOR FROM THu PLANET TO Tht EARTH, 

CX,CY,CZ - UNIT VECTOR FROM THE PLANET TO CANOPUS. 

REAL JC 

RD=57.295779513C823 

CALL ttARTril JD ,XHE ,YFE ,ZHt ,OXri£,DYHE,OZHE > 

CALL PR£Ci=S (243 3282 C6034059 2 , .60342839 ,-. 795 1 3092 , JD ,C Xc , CYE , C 

IZE) 
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APPENDIX D - Concluded 


c 

IF(1B00Y.£Q.4) GO TC 2 

2 CALL £MARS( JD,XhP,YHP»ZHP»DXHP.OYHP,DZHP) 

C 

3 XHPL=XH£-XHP 
YHPL=Yh£-YhP 
ZhP£=ZH£-ZHP 

RSi£=S0Rr(XH£*^2 4YHE^*2 + ZhE*^2) 

RSP=SQRTi XhP«^2+YHP**2+ZHP**2) 
RPL=S 0 RT(XHPE*=^«^ 2 +YHPE »*2 + ZHPe*<‘ 2 ) 

SEX=XhL/RSE 

ScY=Ynt/RSr: 

SL Z = ZHL/RSfc 

SPX=XHP/RSP 

SPY=YHP/RSP 

SPZ=ZHP/RSP 

PLX=XI-Pc/KPL 

PLY=YHPc/RPt 

PEZ=ZhP;:/RP£ 

CALL LAFLNGC SlX,SlY tSEZ tEHLATt EHLONG) 

CALL LATLNG(SPX,SPY tSPZt PHLAT »PHLONG) 

CALL DCT( SEXiSEYf SEZfSPX,SPY,SPZ,ESP) 

CALL DCTlSfcX,S£Y,SEZf PEX,PEYfPEZf SEP) 

CALL DCT(SPX,SPYtSP2»-PEX,-PEY»-PEZfSP6) 

CALL RECeC( JD»-SPX,-SPY,-SPZ,SX£tSYE,SZE) 

CALL RuCLO( JDfPEX,PEY,PLZ*EXEf EYE»EZE) 

C 

IF ( IBGCY.L0.4) GQ TC 5 

5 CALL RECMEO( JOtSXfc ,SYE,SZE»SX,SY,SZ,DfcCSfRAS) 
CALL Rc:CMcO( JCtLXLf EYEfEZEf EXf £Y,EZf OECEtRAE) 
CALL HLQMcQi JCtCXLf CY£,CZ£,CX,CY,CZtDECCfRAC) 

C 

6 XPS^SX’J'RSP 
YPS^SY^^'RSP 
ZPS=SZ^RSP 
XPL=EX«RP£ 

YPE=LY*RPE 

ZPa=:cZ^R.PE 

C 

800 RL TURN 
cKD 
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